Transcriptome Profiling Based on Different Time Points After Hatching Provides a Core Set of Gene Resource for Understanding Larval Immune Response Mechanisms Against Vibrio anguillarum Infection in Amphioctopus fangsiao

Immune defense systems are indispensable for living organisms. Within an immune network, problems with any given link can impact the normal life activities of an organism. Amphioctopus fangsiao is a cephalopod that exists widely throughout the world’s oceans. Because of its nervous system and locomotive organs, it has become increasingly studied in recent years. Vibrio anguillarum is one of the most common pathogenic bacteria in aquaculture organisms. It is highly infectious and can infect almost all aquaculture organisms. V. anguillarum infection can cause many adverse biological phenomena, including tissue bleeding. Study the immune response after V. anguillarum infection would help us to understand the molecular mechanisms of immune response in aquaculture organisms. In this research, we infected the primary incubation A. fangsiao with V. anguillarum for 24 h. We analyzed gene expression in A. fangsiao larvae via transcriptome profiles at 0, 4, 12, and 24 h after hatching, and 1,385, 734, and 6,109 differentially expressed genes (DEGs) were identified at these three time points. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were used to identify immune-related DEGs. Protein–protein interaction networks were constructed to examine interactions between immune-related genes. Twenty hub genes involved in multiple KEGG signaling pathways or with multiple protein–protein interaction relationships were identified, and their differential expression verified by quantitative RT-PCR. We first studied V. anguillarum infection of A. fangsiao larvae by means of protein–protein interaction networks. The results provide valuable genetic resources for understanding immunity in molluscan larvae. These data serve as a theoretical basis for the artificial breeding of A. fangsiao.

Immune defense systems are indispensable for living organisms. Within an immune network, problems with any given link can impact the normal life activities of an organism. Amphioctopus fangsiao is a cephalopod that exists widely throughout the world's oceans. Because of its nervous system and locomotive organs, it has become increasingly studied in recent years. Vibrio anguillarum is one of the most common pathogenic bacteria in aquaculture organisms. It is highly infectious and can infect almost all aquaculture organisms. V. anguillarum infection can cause many adverse biological phenomena, including tissue bleeding. Study the immune response after V. anguillarum infection would help us to understand the molecular mechanisms of immune response in aquaculture organisms. In this research, we infected the primary incubation A. fangsiao with V. anguillarum for 24 h. We analyzed gene expression in A. fangsiao larvae via transcriptome profiles at 0, 4, 12, and 24 h after hatching, and 1,385, 734, and 6,109 differentially expressed genes (DEGs) were identified at these three time points. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were used to identify immune-related DEGs. Protein-protein interaction networks were constructed to examine interactions between immune-related genes. Twenty hub genes involved in multiple KEGG signaling pathways or with multiple protein-protein interaction relationships were identified, and their differential expression verified by quantitative RT-PCR. We first studied V. anguillarum infection of A. fangsiao larvae by means of protein-protein interaction networks. The results provide valuable genetic resources for understanding immunity in molluscan larvae. These data serve as a theoretical basis for the artificial breeding of A. fangsiao.

INTRODUCTION
As a new research direction, mollusks have received extensive attention recently because of their complex immune responses and molecular mechanisms (Sokolova, 2009;Wang et al., 2019). organism growth is inseparable from immune response, and loss of immunity leads to disease and even death, suggesting that immune response is a core factor of the growth of organisms (Rowley and Powell, 2007;Barcia and Ramos-Martínez, 2011;Wang et al., 2019). The larval is a relatively vulnerable stage in life (Ginger et al., 2013;Kaplan et al., 2013). Growth quality among larvae directly affects the quality of life of the adult, impacting characteristics including vitality and physique. The immune functions of larvae are important to the growth process. It is well known that organisms have a higher ability to resist invasion by external pathogens when the numbers and activity of their immune cells are high. Therefore, immune function can effectively protect larvae of organisms, maintaining their normal growth (Wang et al., 2013;Matozzo, 2016).
As a marine mollusk, Amphioctopus fangsiao is widely distributed along the coasts of the Pacific, the Yellow Sea, and the Bohai Sea. Due to the difficulties associated with catching wild A. fangsiao, numbers of A. fangsiao bred in captivity have gradually risen in recent years. However, whether A. fangsiao is wild or bred the larvae are vulnerable to a variety of pathogens including bacteria and viruses (Budelmann, 1994;Wei et al., 2015Wei et al., , 2018. Vibrio anguillarum is a grampositive bacterium that can cause many popular diseases in aquaculture animals. In recent years, multiple studies have confirmed V. anguillarum infections in many fish species, including Salmo salar and Paralichthys olivaceus (Hoel et al., 1997;Li et al., 2015;Bao et al., 2019). V. anguillarum can infect bivalve and gastropod such as Chlamys nobilis and Haliotis gigantea (Cong et al., 2008;Iehata et al., 2009;Zhang et al., 2016). When aquaculture organisms are infected by V. anguillarum, they display red gills, muscle bleeding, mucosal tissue decay, and other symptoms, ultimately leading to significant losses to the aquaculture industry. V. anguillarum infection requires certain conditions, and the incidence rate of V. anguillarum infection tends to be higher in poor environments or among injured or stimulated organisms (Vázquez et al., 2006;Zhou et al., 2010;Nie et al., 2017). Although there are widely studies on the infection of V. anguillarum in fish, bivalve, and gastropod, there are still few studies have been conducted on the infection of mollusks, especially cephalopods. The pathogenic mechanisms of V. anguillarum in cephalopods are still unclear, and require further study.
Transcriptome sequencing technology has recently become an important method for studying differences between individuals of the same species. To aid in assessing biological immune response mechanisms, transcriptome sequencing can analyze differences in gene expression between different individuals or tissues (Zhao et al., 2012;Silva et al., 2013;Zhang et al., 2018). It can help identify important immune genes and further study biological immune response mechanisms. Previous researches have carried out transcriptome analyses of aquatic mollusks such as Pinctada martensii (Zhao et al., 2012), Biomphalaria glabrata (Silva et al., 2013), Haliotis diversicolor (Zhang et al., 2018), Octopus vulgaris (Castellanos-Martínez et al., 2014), and Euprymna tasmanica (Salazar et al., 2015;Schultz and Adema, 2017). The transcriptome analyses of A. fangsiao is still lacking, which need a lot of experiments to explore.
In this research, we infected the primary incubation A. fangsiao larvae with V. anguillarum for 24 h. We then took A. fangsiao larvae at 0, 4, 12, and 24 h after hatching for transcriptome sequencing and bioinformatics analyses, including gene function annotation, differentially expressed genes (DEGs) analysis, Gene Ontology (GO) functional enrichment analysis, Kyoto Encyclopedia of Genes and Genomes (KEGG) functional enrichment analysis, and immune-related protein-protein interaction network studies. We identified and verified 20 differentially expressed hub genes using quantitative RT-PCR. These results revealed the immunity of A. fangsiao larvae after infection with V. anguillarum, providing a valuable resource for understanding immune response mechanisms in mollusks, and laying a foundation for exploring the immune response mechanisms of other mollusks.

Ethics Statement
Amphioctopus fangsiao samples were obtained from a commercial hatchery. This research was conducted in accordance with the protocols of the Institutional Animal Care and Use Committee of the Ludong University (protocol number LDU-IRB20210308NXY)

Sample Collection and RNA Preparation
The wild parent used in these experiments were obtained from the Rizhao Sea area. After spawning, the eggs were protected and incubated by female parent A. fangsiao. During the incubation process, the seawater temperature was maintained at 19-20.8 • C, and the incubation lasted approximately 29 days. After hatching, larvae were temporarily cultured in floating seawater containing 1 × 10 7 CFU/mL V. anguillarum for 24 h. Primary incubation larvae were collected and stored for 4, 12, and 24 h in liquid nitrogen until RNA was extracted using the TRIzol method.
Nine larvae from primary incubation were randomly selected from each time point for RNA extraction: A. fangsiao larvae infected for 0 h (Oo-C), A. fangsiao larvae infected for 4 h (Oo-4h), A. fangsiao larvae infected for 12 h (Oo-12h), and A. fangsiao larvae infected for 24 h (Oo-24h). At each time point, equal molar masses of RNA from any three of nine larvae were pooled into one replicate to provide a template for the construction of the transcriptome library; any three of the remaining six larvae were pooled into the second replicate; and the final three larvae were pooled into the third replicate. Remaining RNA was retained for subsequent quantitative RT-PCR verification.

Library Construction and Illumina Sequencing
Library construction was carried out according to a method described by Li et al. (2017). Finally, we sequenced each sample using an Illumina HiSeq 4000 platform.

Gene Expression Level Analysis and Functional Annotation
A reference sequence was obtained by splicing Trinity with RSEM, and clean reads were mapped to this reference sequence. We then obtained the number of read counts for each sample that mapped to each gene. The read count was transformed into fragments per kilobase of transcript per million mapped reads (FPKM) to analyze the gene expression and abundance, and expression levels of samples were positively correlated with FPKM. The unigenes were annotated by searching the sequences against the NR, NT, GO, SwissProt, and KOG databases using BLASTX with a cut-off of E value ≤ 1e-5. Meanwhile, we annotated unigenes into Pfam database with Hmmer 3.0 package (E value ≤ 0.01).

Differentially Expressed Genes Screening and Analysis
Differentially expressed genes were screened out using the DESeq2 package for R as a negative binomial distribution model. First of all, data were imported for building the dds model, and then the DESeq function was used to estimate the dispersion of the samples. After that, the difference in gene expressions was analyzed by this package. Finally, DEGs with q-value ≤ 0.05 and |log2 fold change| ≥ 1 were screened out (Love et al., 2014;Ho et al., 2017). Significant DEGs in each time point were screened by GO functional enrichment analysis to obtain GO terms and distribution of DEGs. Statistical analyses of DEGs were conducted. To further understand functions of DEGs, KEGG signaling pathway analysis was performed for up-and down-regulated genes, and we annotated immune pathways involving DEGs. Finally, we screened the signaling pathways in which DEGs were significantly enriched, and species and numbers of DEGs in these pathways were counted.

Trend Analysis
Through trend analysis, we understood the change trend of gene expression over time. STEM was used to analyze the expressing trend of DEGs, and cluster trends in gene expression quantity. The parameters were set as follows: the Maximum Unit Change in model profiles between time points is 1; maximum output profiles number is 20 (similar profiles will be merged); minimum ratio of fold change of DEGs is no less than 2.0; the criterion for screening trends with significant differences is p-value ≤ 0.05 (Ernst and Bar-Joseph, 2006). Trend analysis showed the expressing trend of genes in the process of larval growth after infection, and predicted the future trend, which can help us explore the relationship between immunity and larval growth of A. fangsiao after infection.

Functional Protein Association Network Construction
We used STRING v11.0 with default parameters to construct protein-protein interaction networks to further explore gene relationships in immune regulation pathways (Damian et al., 2018).

Quantitative RT-PCR Validation
We verified the accuracy of our RNA-Seq results using quantitative RT-PCR to validate 20 selected genes (Supplementary Table 1). Three biological replicates of each time point were used in the experiment. Primer Premier 5.0 was used to design gene-specific primers. Table 1 lists the 20 selected genes and their corresponding primer sequences. We evaluated the stability of GAPDH, β-actin, and 18S genes in different tissues and embryo development stages of A. fangsiao. The expression level of A. fangsiao β-actin tended to be stable, and it was used as an endogenous control in this experiment. Quantitative RT-PCR was performed according to a method described by Li et al. (2019).

Sequencing Results and Quality Assessment
Samples at four time points were sequenced by RNA-Seq. The overview of sequencing results is provided in Table 2. Raw sequencing reads were submitted to Sequence Read Archive in NCBI; the SRA accession numbers were SRR15204591, SRR15204592, SRR15204593, SRR15204594, SRR15204595, SRR15204596, SRR15204597, SRR15204598, SRR15204599, SRR15204600, SRR15204601, and SRR15204602 2 .

Differential Expression Analysis
Through differential expression analysis, we identified 1,385, 734, and 6,109 DEGs at 4, 12, and 24 h after infecting compared with primary incubation larvae. Among these, 480 DEGs were up-regulated, and 905 DEGs were down-regulated at 4 h after infection; 338 DEGs were up-regulated, and 396 DEGs were down-regulated at 12 h after infection; 1,685 DEGs were upregulated, and 4,424 DEGs were down-regulated at 24 h after infection (Figure 1 and Supplementary Tables 2-4). A Venn diagram shows all DEGs after infecting. Because DEGs at each of the three time points might be key genes affecting immunity in A. fangsiao larvae, their union set (7,019 DEGs) was selected for subsequent analysis (Figure 2 and Supplementary Table 5). The heatmap (Figure 3 and Supplementary Table 6) intuitively shows the clustering distribution of these DEGs. Differentially Expressed Genes Expression Trend STEM was used to analyze trends in 7,019 DEGs. The results showed that DEGs were divided into 20 expression trends, involving six significantly enriched ones (Figure 4). Their p-value were less than 0.05, and 4,940 DEGs were differentially expressed in these trends. Among them, profile 0 ( Figure 4D) is the most significant trend, and profile 9 ( Figure 4C) has the most genes expressed in this trend.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes Enrichment Analysis of Differentially Expressed Genes
We identified the top ten level-3 terms of the three categories (biological process: 352 level-3 subclasses, molecular function:   Table 7) by GO functional enrichment analysis. For instance, the signal transduction term and the cell differentiation term in biological process were closely related to immunity. Similarly, KEGG enriched DEGs in different pathways guided us to further understand the functions of these genes. Of 7,019 DEGs, 1,548 DEGs were enriched in 263 level-2 KEGG class pathways (Figure 6 and Supplementary Table 8). Table 3 showed the 20 immune-related pathways that were significantly enriched including biosynthesis of antibiotics signaling pathway, Wnt signaling pathway, and PI3K-Akt signaling pathway. Because some DEGs were involved in the regulation of multiple signal pathways, we finally obtained 208 genes enriched in 20 immune-related pathways for follow-up analyses.

Construction of Immune-Related Protein Interaction Networks
Proteins are important components of all cells and tissues in living organisms. All important parts of an organism require functional proteins, which are the main mediators of the activities of life. Construction of protein-protein interaction networks can help reveal the key genes involved in immune processes. We constructed protein-protein interaction networks using the protein sequences encoded by 208 genes in 20 significantly enriched immune-related signaling pathways. Figure 7 shows these networks. Among this figure, the CASP gene had the highest number of protein-protein interactions, and the other 5 genes including AXN1, CREBBP, MTOR, PTK2, and WNT5A were interacted with more than 20 proteins. Relevant parameter information of them was listed in Table 4. There were 148 nodes and 622 edges in the network, and each node was interacted with an average of 8.41 nodes. The clustering coefficient was 0.466, and the p-value of this network was ≤ 1.0e-16.

Acquisition and Validation of Key Differentially Expressed Genes Related to Immune Processes
This research focused on interactions between key immunerelated genes. Twenty key DEGs (Table 5) involved in multiple signaling pathways or multiple interaction relationships were obtained from KEGG pathway analysis and protein-protein interaction networks. These were then used to explore interaction mechanisms between them. The 20 key genes were divided into five categories as follows: Collagen family, Wnt signaling pathway, biosynthesis of antibiotics signaling pathway, tumor necrosis factor (TNF) signaling pathway, and other important immune-related genes. To verify the accuracy of these key genes, we used quantitative RT-PCR to quantitatively assess the relative expression of 20 immune-related genes at each time point identified as DEGs by RNA-Seq. Quantitative RT-PCR results showed that all DEGs measured were single products. A comparison of gene expression profiles by RNA-Seq and quantitative RT-PCR (Figure 8) indicated that the quantitative RT-PCR results correlated significantly with the RNA-Seq results, and the two methods presented the same trend pattern.

The Purpose and Significance of This Research
Amphioctopus fangsiao has become increasingly popular because of its rich nutritional value and delicious meat. As an economically important species employed in mariculture, the artificial breeding of A. fangsiao has been well concerned (Iehata et al., 2009;Zhang et al., 2016;Li et al., 2019). The larval growth stage is very important for A. fangsiao, which determines whether the adult can grow healthily. The larvae of aquaculture organisms are easily infected by a variety of pathogens, which can seriously hinder aquaculture development (Barcia and Ramos-Martínez, 2011;Wang et al., 2019). Research into biological immune responses to bacterial infection will help breeders to prevent and mitigate bacterial infections. We infected primary incubation larvae for 24 h with V. anguillarum, and identified the changes of the biological immune response mechanisms. We identified 7,019 DEGs and believed that these genes have an important relationship with the immune response mechanisms of A. fangsiao larvae. A heatmap showed that there were significant differences in immunity with increase in A. fangsiao larvae infection time. The results of trend analysis showed that the six of the 20 trends were significantly enriched, and half of DEGs were expressed in these trends, indicating significant differences in the immunity after larval infection. KEGG enrichment analysis obtained 20 immune-related pathways enriched in DEGs, which were interacted with each other and shared biological immune functions. Finally, we used 208 genes in these pathways to construct protein-protein interaction networks to understand immune response mechanisms of A. fangsiao larvae.

Trend Analysis of Differentially Expressed Genes
The results of the trend analysis in Figure 4 indicated that profile 9 trend had the largest number of DEGs, and the significance of profile 0 trend was the highest. Among them, the gene expression trends showed continuously down-regulation and from stable to down, respectively, which were inconsistent with the expectation of our experiment. Therefore, we conducted GO and KEGG enrichment analyses of the genes involved to further explore the causes. In our analysis, we found that these genes were mainly involved in metabolic processes, including protein digestion and absorption signaling pathway, insulin signaling pathway, and cell adhesion term (Curtis et al., 1978;Lee and Pilch, 1994;Schepis et al., 2012;Brezas and Hardy, 2020). These results showed that the metabolic capacity of A. fangsiao larvae decreased gradually in the course of infection. Meanwhile, 1,291 DEGs in profile 17 and profile 19 were upregulated. We used DEGs from these two profiles for GO and KEGG enrichment analyses, and identified a number of immunerelated terms and pathways. For instance, positive regulation of JNK cascade term, biosynthesis of antibiotics signaling pathway, and NF-kappa B signaling pathway were significantly enriched. The JNK cascade is a key step in cell apoptosis mediated by the NF-kappa B signaling pathway, which is involved in the apoptosis process of diseased cells (Bubici et al., 2004;Lanna et al., 2017). NF-κB is involved in the regulation of immunity, lymphocyte development, and tumorigenesis (Bubici et al., 2004;Sun, 2017). Antibiotics can effectively resist the invasion of bacteria, and their synthesis can improve biological resistance (Bulfon et al., 2020). In conclusion, although the DEGs for the two most critical trends were down-regulated, they were highly correlated with metabolism. The results indicated that the metabolism of larvae was affected and had a downward trend after infection. The other DEGs in two up-regulated trends were more associated with immunity, suggesting that the immune response of A. fangsiao larvae increased gradually with the prolongation of infection. Above results revealed significant changes in the responses of A. fangsiao after infection. The complex mechanisms of A. fangsiao after infection need to be further studied and discussed in the future.

Enrichment Analysis of Immune-Related Kyoto Encyclopedia of Genes and Genomes Signaling Pathways and Gene Ontology Terms
Many immune-related terms and pathways were yielded through GO and KEGG enrichment analyses, including signal transduction term, cell differentiation term, biosynthesis of antibiotics signaling pathway, TNF signaling pathway, and other key terms and pathways. These terms and pathways were vital in biological immunity, indicating that there were a large number of immune-related genes in the larvae infected with V. anguillarum, and thus resulted in complex immune responses in the larvae. A comprehensive analysis of the enriched KEGG pathways and GO terms contributes to investigate the molecular mechanisms and immune systems changes of A. fangsiao larvae after infection with V. anguillarum, thereby aiding our artificial breeding work.

Speculation of Hub Genes
Proteins are the basis of biological activities, and their interactions maintain the growth and development of organisms. We used 208 key DEGs identified in immune-related pathways to construct protein-protein interaction networks. The results suggested that interactions between these proteins was greater than expected relative to a randomly selected a group of proteins of similar size from the genome. This conclusion proves that above proteins cooperate in carrying out coordinated functions. We thus suggested that nodes with more edges were hub proteins in immune response. The genes corresponding to these hub proteins were identified as hub genes for further studies and verification.

Functional Analysis of Kyoto Encyclopedia of Genes and Genomes Signaling Pathways and Hub Genes
Transcriptome profiling was analyzed to compare the changes of immune response mechanisms in A. fangsiao larvae infected with V. anguillarum within 24 h. This helps us to further understand the A. fangsiao larval immune response to V. anguillarum infection with the growth of them. Finally, we speculatively identified 20 hub genes involved in multiple protein-protein interaction relationships or KEGG signaling pathways, and these hub genes, KEGG signaling pathways, and immune differences at each time point were investigated.

Collagen Family
Collagen is one of the most abundant proteins in the organisms. It is an important part of extracellular matrix, and plays a regulatory role in tissue structure. Its subunits directly affect cell phenotypes by mediating cell differentiation and apoptosis via relevant signaling pathways (Brondijk et al., 2010;Jürgensen et al., 2020). Collagen also acts as a major part in the immune system. When cells are stimulated by bacteria, collagen receptors are responsive to change the extracellular environment, thus indirectly regulating the activities of immune cells to resist the invasion of bacteria (Kresina et al., 1984;Aragona et al., 1998;Brondijk et al., 2010). In this study, five key genes involved in the immune response of collagen family (COL1A1, CoL1A2, COL2A1, COL4A1, and COL4A4) were identified. Among them, COL1A1 and COL2A1 presented the most obvious effect on enhancing the biological resistance to external invasion through the protection of bones, kidney, skin, and other organs. The lack of collagens in the organisms can lead to osteogenic and skin damages, which can reduce the immunity of organisms indirectly (Kuppevelt et al., 1995;Takai et al., 2001;Makino et al., 2010). In our study, the five genes in collagen family were enriched significantly, indicating an obvious bacteriological stimulation effect. Collagen family may be an important part to regulate the biological immune functions. In the above five genes, COL1A1 and COL2A1 were significantly up-regulated with the prolongation of infection time. This result indicated that A. fangsiao larvae had a strong immunity to effectively resist the invasion of pathogens. The other three genes were continuously downregulated within 24 h of infection, the specific reasons need additional explorations.

Wnt Signaling Pathway
Wnt is an important immune regulatory protein widely existing in various organisms. It participates in the immune response of multiple organs and acts a key role in maintaining tissue homeostasis (Famili et al., 2016). The Wnt signaling pathway can be divided into two subtypes: typical Wnt signaling pathway and atypical Wnt signaling pathway. The interaction between these two pathways regulates the growth, proliferation, differentiation, and apoptosis of immune cells, which can affect multiple immune responses of organisms (Naskar et al., 2014;Yu et al., 2014;Chae and Bothwell, 2018). Moreover, the Wnt signaling pathway is responsible for maintaining homeostasis through repairing damaged tissues and regulating the inflammatory response, which plays a role in the resistance and treatment of some immune diseases (Staal et al., 2008;Pai et al., 2017;Ruan et al., 2021). In this study, the Wnt signaling pathway was significantly enriched, and multiple DEGs including TBL1XR1, ROCK2, and CSNK1A1 were enriched in this pathway. They are the key factors in the immune system, and play a significant role in immune regulation during larval growth. Meanwhile, these genes are functional in the resistance of larvae against the invasion of pathogens and promote the normal growth of organisms (Puigvert et al., 2013;Nikola et al., 2017;Venturutti et al., 2020). Their up-regulation after larval infection suggested the improved immune ability of A. fangsiao larvae after being infected with bacteria.

Biosynthesis of Antibiotics Signaling Pathway
Antibiotics are a class of metabolites produced in the life of higher animals, plants, and fungi, which have the ability to resist the invasion of pathogens. They can affect the growth of organisms by regulating cell growth and apoptosis (Pomorska-Mól and Pejsak, 2012;Herman, 2020). Biosynthesis of antibiotics can improve the ability of organisms to resist bacterial infection and increase the resistance of larvae. When organisms are invaded by FIGURE 7 | Immune-related protein-protein interaction networks. Network nodes stand for proteins. The legend represents the relationships between nodes.
Frontiers in Marine Science | www.frontiersin.org  external bacteria, antibiotics participate in the immune response, which can help organisms resist invasion, and increase their immunity. Recent studies have shown that antibiotics can treat some immune diseases and prevent immune system disorders (Stevens, 1996;Pomorska-Mól and Pejsak, 2012;Yang et al., 2017). In our study, biosynthesis of antibiotics signaling pathway was significantly enriched, suggesting that antibiotics act as an important part in the immune response of A. fangsiao larvae. Many key immune-related genes in this pathway were up-regulated with the growth of infected A. fangsiao larvae, indicating that antibiotics can guide a variety of immune responses to eliminate invading pathogens.

Tumor Necrosis Factor Signaling Pathway
As a regulatory factor, TNF exists in various organisms. It can send signals to proteins in the immune system, promote the interaction between immune cells, and enhance the immune ability of organisms (Ferreira et al., 1993;Vielhauer and Mayadas, 2007;Croft, 2009;George et al., 2014;Li et al., 2021). In the immune system, TNF also acts as a signal protein to regulate the release of other cytokines. Moreover, TNF plays an important role in promoting T cell proliferation, interferon expression, and accelerating the efficiency of chemokines (Palanki et al., 2000;Vielhauer and Mayadas, 2007;Lougaris et al., 2017;Li et al., 2021). In our research, multiple immune-related genes including NFKB1, ATF4, and MAP2K4 were significantly enriched in this pathway. NFKB1 plays a major part in the regulation of immune response in organisms, which increases the resistance of organisms by regulating expression levels of proinflammatory factors, and inflammation-related proteins in immune cells. The other two genes are also useful in regulating the immune response of organisms (Sethi et al., 2007;Vielhauer and Mayadas, 2007;Huang et al., 2019). This result suggested that immune genes responded quickly after infection, which increasing the immune activity of A. fangsiao larvae, and effectively eliminated the invading pathogens.

Analysis of Three Important Hub Genes
Among all hub genes, CASP3, CREBBP, and PTK2 were enriched in more than two immune-related signaling pathways, and possessed higher numbers of protein-protein interactions. CASP3 is a regulator of cell apoptosis. It can maintain the dynamic balance of cell number and the healthy growth of organisms by regulating the apoptosis of pathological cells and aging cells. The up-regulation of CASP3 is helpful to clear the diseased cells, and improve the survival rate and vitality of organism larvae (Lowe and Lin, 2000;Suzuki et al., 2020).
As an acetyltransferase, CREBBP exists widely in lymphocytes of various organisms. It is a key gene regulating immune response, resisting pathogen invasion, and promoting immune cell expression. The low expression or loss of function of CREBBP leads to the disorder of enhancer network and B cell signal, and causes the disorder of expression of other immune cells (Jiang et al., 2016;Zhang et al., 2017;Jia et al., 2018). PTK2 is a protein kinase widely existing in various tissues, and organs of organisms, which acts as an extremely important disease regulator. At the same time, PTK2 is involved in the regulation of cell adhesion and movement, which is conducive to substance transport, and the expression levels of immune genes Yi et al., 2021). In this study, CREBBP was up-regulated within the first 12 h of infection, and slightly down-regulated within the last 12 h, indicating an important role in the immune regulation of A. fangsiao larvae after infection. PTK2 was significantly upregulated within 24 h after infection and down-regulated within 4 to 12 h. The causes of PTK2 down-regulation need to be explained by subsequent studies. The above two genes had obvious effects on larval immunity, which is beneficial for the larvae to resist the invasion of pathogens. However, with the prolongation of infection time, CASP3 was down-regulated. This result is not consistent with our study, which needs additional studies to explain.

Other Hub Genes and Kyoto Encyclopedia of Genes and Genomes Signaling Pathways
In addition to the above hub genes and KEGG signaling pathways, our results identified other key DEGs and pathways involved in A. fangsiao larval immunity, including pathways in cancer signaling pathway, NF-kappa B signaling pathway, the PLCG1 gene, the ITGA4 gene, and the MTOR gene. The above two signaling pathways and three genes have been reported to be closely related to the immune response in various organisms Thurnher et al., 2013;Darzi et al., 2017;Dieli-Crimi et al., 2018;Hong et al., 2019;Shiseki et al., 2019). The signaling pathways and genes details of the immune response mechanisms of A. fangsiao larvae infected with V. anguillarum require additional studies to further elucidate mechanistic details.

CONCLUSION
We performed transcriptome analyses of gene expression in primary incubation A. fangsiao larvae within 24 h of V. anguillarum infection, and established a protein-protein interaction network. Twenty hub DEGs enriched in multiple KEGG signaling pathways or with multiple protein-protein interaction relationships were identified. In our research, protein-protein interaction networks were first used to explore differences on A. fangsiao larval immunity. Our results provide valuable genetic resources for understanding the immunity of invertebrate larvae. Meanwhile, the data serve as a foundation for additional research into immune processes in invertebrates.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are publicly available in NCBI using accession numbers SRR15204591 -SRR15204602 at the following link: https://www.ncbi.nlm.nih. gov/Traces/study/?acc=SAMN20309867&o=acc_s%3Aa.