Comprehensive Relationship Analysis of the Long Noncoding RNAs (lncRNAs) and the Target mRNAs in Response to the Infection of Edwardsiella anguillarum in European eel (Anguilla anguilla) Inoculated with Freund’s Adjuvant

Freund’s complete adjuvant (FCA) and incomplete adjuvant (FIA), generally applied in subunit fishery vaccine, have not been explored on the molecular mechanism of the non-specific immune enhancement. As long noncoding RNAs (lncRNAs) play vital regulating roles in various biological activities, in this study, we examined the genome-wide expression of transcripts in the liver of European eel (Anguilla anguilla, Aa) inoculated with FCA and FIA (FCIA) to elucidate the regulators of lncRNAs in the process of Edwardsiella anguillarum (Ea) infection and Aa anti-Ea infection using strand-specific RNA-seq. After eels were challenged by Ea at 28 days post the first inoculation (dpi), compared to the control uninfected eels (Li group), the control infected eels (Con_Li group) showed severe bleeding, hepatocyte atrophy, and thrombi formed in the hepatic vessels of the liver, although eels inoculated with FCIA (FCIA_Li group) also formed slight thrombi in the hepatic vessels. Compared to the FCIA_Li group, there was about 10 times colony-forming unit (cfu) in the Con_Li group per 100 μg liver tissue, and the relative percent survival (RPS) of eels was 50% in FCIA_Li vs Con_Li. Using high-throughput transcriptomics, differential expressed genes (DEGs) and transcripts were identified and the results were verified using fluorescence real-time polymerase chain reaction (qRT-PCR). Interactions between the differential expressed lncRNAs (DE-lncRNAs) and the target DEGs were explored using Cytoscape according to their co-expression and co-location relationship. We found 13,499 lncRNAs (10,176 annotated and 3423 novel lncRNAs) between 3 comparisons of Con_Li vs Li, FCIA_Li vs Li, and FCIA_Li vs Con_Li, of which 111, 110, and 129 DE-lncRNAs were ascertained. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of DEGs targeted by DE-lncRNAs revealed these DEGs mainly involved in single-organism cellular process in BP, membrane in CC and binding in MF, and KEGG pathways showed that the target DEGs in co-expression and co-location enriched in cell adhesion molecules. Finally, 118 DE-lncRNAs target 1161 DEGs were involved in an interaction network of 8474 co-expression and 333 co-location-related links, of which 16 DE-lncRNAs play vital roles in anti-Ea infection. Taken together, the interaction networks revealed that DE-lncRNAs underlies the process of Ea infection and Aa anti-Ea infection.


Introduction
Edwardsiella anguillarum (Ea) is a common bacterial pathogen mainly infected cultivated eel (Reichley et al. 2017;Guo et al. 2019). At present, the research on Edwardsiella focuses on finding a protein with universal antigenicity to different serotypes. The results of some researchers found a universal bivalent expressed outer membrane protein (OMP) generally had good immunogenicity (Guo et al. 2013;Zhang et al. 2016). A study from earlier stage revealed an expressed outer membrane protein OmpS 2 showed good immunogenicity and can be applied as a kit for E. tarda detection (Kumar et al. 2007). Simultaneously, an outer membrane protein A (OmpA) from Ea was also evaluated as an antigen with good immunogenicity in Japanese Flounder Paralichthys olivaceus (Tang et al. 2010) and Anguilla anguilla (Duan et al. 2019). Additionally, examination of the transcription profile of immune-related genes showed that vaccination with OMP (E. tarda 49) induced the expression of a broader spectrum of genes that are likely to be involved in humoral and innate cellular immunity (Jiao et al. 2010).
The above-mentioned OMP vaccine were emulsified with Freund's complete adjuvant (FCA) or incomplete adjuvant (FIA) to enhance the protective effect, and the control adjuvant group designed by researchers also produced protective effect after the challenge of different bacterial pathogens. The relative percentage survival (RPS) for FIA-adjuvanted vaccine (77.8%) was higher than the vaccine without adjuvant (40.7%) after FIA incorporated with inactivated Streptococcus agalactiae for Nile tilapia fingerlings (Ramos-Espinoza et al. 2020) as well as the FCA adjuvant induced non-specific immunity in cultured fish and increased 5.3 to 450 times of LD 50 on Aeromonas salmonicida, A. hydrophila, and Vibrio ordalii challenge in coho salmon Oncorhynchus kisutch, respectively (Olivier et al. 1985).
However, there is no study focused on FCA and FIA inoculated sequentially and no molecular mechanism of the non-specific immune protection activated by the FCA and FIA (He et al. 2020a, b). Recently, it was found that the Ea mainly invaded the liver of European eels post 48 h of the infection (Zhai et al. 2021), and two times of immunization of OmpA emulsified with FCA and FIA respectively showing long non coding RNA (lncRNA) played crucial role in regulating many immune-related proteins (He et al. 2021). LncRNA, a kind of RNA molecule with a length of more than 200 nt, have been ignored as useless transcripts for a long time (Yousefi et al. 2020). In recent years, more and more studies have proved that lncRNA is widely existed in eukaryotes, and plays an important role in all kinds of animals (Robinson et al. 2020;Valenzuela-Muñoz et al. 2021).
In fish, integrated analysis of lncRNA and mRNA expression showed anti-infection regulation of lncRNAs post the infection of E. tarda and Ea (Xiu et al. 2021;Zhao et al. 2021;He et al. 2021). However, there is no study on the anti-infection effect of lncRNAs activated by the bacterial pathogen after animals inoculated with adjuvant. Since the adjuvants of FCA and FIA have been widely used in fishery vaccine, it is of great significance to explore its molecular immune mechanism. The aim of this study is to explore the relationship between the differential expressed lncRNAs (DE-lncRNAs) and the target DE-mRNAs activated by the Ea infection pre and post the inoculation of FCA and FIA in European eels (Anguilla anguilla, Aa). Results of this study will reveal the molecular response of Freund's adjuvant to bacterial infection in fishery vaccine.

Bacterial Strain and Fish
The pathogen of Edwardsiella anguillarum (Ea B79) was isolated from the liver of diseased eels and identified by the Biolog GeneIII system and sequencing analysis of OmpA genes (Duan et al. 2019). One hundred and twenty European eels (40 ± 10 g) were obtained from QingLiu county, Fujian province, and randomly assigned into 3 groups. Eels were placed in 500-l tanks under controlled conditions with aerated fresh water and kept at 28 ℃ for 14 days prior to the inoculation. During the experiment, eels were fed commercial pellets (Tianma Feed Co., Ltd, Fuqing, China) once a day. Feces of fish and feed remnants were removed after 2 to 3 h of feeding.

Inoculation, Experimental Challenge, and Liver Sampling
Culture of Ea was prepared using Tryptone soya broth (TSB), and the bacterial cells were concentrated to 3.0 × 10 7 cfu ml −1 . The Freund's complete adjuvant (FCA) and incomplete adjuvant (FIA) were mixed with equal volume of PBS (0.01 M, pH = 7.4). After the injection of 0.2 ml PBS per fish at 0 day and 14 days, 80 European eels were divided equal into 2 groups of PBS injected (Li group) and Ea infected (Con_Li group) fish which received 0.2 ml/fish of PBS and 3.0 × 10 7 cfu/ml Ea respectively at 28 day post the first injection (dpi). Simultaneously, 40 eels were inoculated with 0.2 ml/fish of FCA and FIA at 0 d and 14 d respectively (FCIA_Li group), and all eels were challenged by the Ea at 28 days post the FCA inoculation. At 24 and 48 h post the challenge (hpc) respectively, 3 eels in each group were anesthetized by immersion in water containing 100 ppm Eugenol for 10 min and then all livers were collected. The remaining eels in two infected groups were used to ascertain the relative percent of survival (RPS) in 14 days post the challenge (dpc) as the method of He et al. (2020a, b).

Bacterial Load, Pathological Observation, and Total RNA Extraction
After 18 livers collected in the method of 2.2 were carefully examined, 9 livers (3 livers/group) collected at 24 hpc were used to detect the colony quantity of bacterial load in the liver (cfu/mg). Another 9 livers (3 livers/group) collected at 48 hpc were cut into 3 pieces. One piece of each liver was placed in 10% medium formalin and processed for paraffin embedding as He et al. (2021), and 4-µm sections were cut from each liver and stained with hematoxylin and eosin (H & E). Images of each H & E-stained liver were acquired by Leica Scan 500 slide scanner equipped with 20 × lens. Another 2 pieces of each liver were collected for bacterial plate colony count and total RNA extraction respectively as the method of Zhai et al. (2021).

Transcriptome of Strand-Specific RNA Sequencing
To obtain the RNA-seq data from 9 liver samples, the flow process of sample detection, database construction, and sequencing in Novogene (Beijing, China) is as follows: total RNA sample detection; rRNA removal; synthesis of double stranded cDNA; cDNA purification; ending cDNA repair and adding A-tail and adaptor; fragment selection, degradation of the second chain of cDNA, PCR enrichment, library quality inspection and illumine 3000 sequencing. To detect the mapping rates, clean reads were aligned with the genome of A. anguilla by the software of Biwtie2 (https:// www. ncbi. nlm. nih. gov/ assem bly/ GCF_ 01334 7855.1).

Genome Alignment and Quantification of Gene and Transcripts
The expression of each gene was quantified by RSEM software. The proportions of reads in exons, introns, and intergenic regions according to the comparison results were calculated. The reads in the intron region may come from the precursor mRNA or the introns retained by alternative splicing events, and reads in the intergenic region may come from lncRNA (Alessandra et al. 2017). The software of Stringtie has specific parameters for different libraries, which can accurately splice transcripts and achieve transcriptional quantification (Pertea et al. 2015). In this study, Stringtie was used to splice the reads which mapped on the genome into transcripts and quantify it (Ghosh and Chan 2016). All transcripts were merged by the Cuffmerge software to remove transcripts length less than 200 nt or with uncertain chain direction. The remaining transcripts were compared with known databases and the gene function annotation of known transcripts were obtained according to the annotation file of A. anguilla (.gff). For unknown transcripts, CPC2, Pfam, and CNCI were used to predict the coding potential. If the transcripts had no coding potential by these software, the candidate transcripts are predicted as novel_lncRNAs (XR_), and if the candidate transcripts are found to have coding potential, it will be predicted as novel_mRNAs (Harrow et al. 2012).

Analysis of Differentially Expressed lncRNA
In this study, up or down functional Gene, mRNA, and lncRNA were classified based on fragments per kilobase per million (FPKM). When compared to 9 samples, the total clusters of DEGs or DE-transcripts found by EdgeR with padj < 0.05 were showed by heat and volcano map (Varet et al. 2016). Using the mainstream hierarchical clustering method, the log10 (fpkm + 1) value was transformed and clustered. The clustering results are represented by Heatmap of Gene, mRNA, and lncRNA. Venn and volcano map, which directly showed the cross and overall distribution genes with significant difference expression between 3 groups, were constructed to find the key lncRNAs in the process of anti-infection. Functional analysis of gene ontology (GO) and Kyoto Encyclopedia of genes and genomes (KEGG) focuses on DEGs with biological functions. GOseq (Young et al. 2010) based on Wallenius non central hyper geometric distribution mathematical model was used to analyze top 20 significant GO enrichments (q value < 0.05) of molecular function (MF), cellular component (CC), and biological process (BP) in DEGs. KOBAS software (version 2.0) was used to analyze the KEGG pathway enrichment of DE-lncRNAs (padj < 0.05) target genes based on hyper geometric test.

RT-PCR Validation
To validate the results of gene expression obtained from RNA-seq, cDNA samples of 23 immune-related genes were selected and verified by qRT-PCR between FCIA_Li and Con_Li. The analysis was performed using the RealUniversal Color PreMix (SYBR Green) (Tiangen, Beijing, China) on a LightCycler® 480 (Roche, USA). The first-strand cDNA of DEGs of Aa was first synthesized as templates for qRT-PCR with specific primers, and β-actin was used as an internal control (Table 1). Preparation of samples and the amplification procedure were referenced as He et al. (2020a, b). Data were analyzed using the Excel and relative expression of DEGs was determined using the 2 −∆∆Ct algorithm, and then compare to the results of RNA_seq.

Interaction Network Between DE-lncRNAs and Target DEGs
Location-related target gene analysis was used by colocation (prediction of cis target gene according to the location relationship of lncRNA-mRNA, the range of screening is within 100 k) and expression-related target gene analysis was used by co-expression (number of samples = 9). Prediction of target gene was according to the expression correlation between DE-lncRNA (padj < 0.05) and mRNA when the correlation coefficient was greater than 0.95. To analyze the profiles of lncRNA-mRNArelated DEGs interaction network and explore the potential roles of key lncRNAs in the infected eel pre versus post Freund's adjuvant inoculation, the correlation of colocation and co-expression analysis of DE-lncRNAs and mRNA-related DEGs were used to analyze the anti-infection related lncRNAs. The Pearson correlation coefficient between lncRNA and mRNA greater than 0.99 was used as target DEGs (pvalue < 0.05) for constructing network map of the lncRNA-mRNA interaction. Cytoscape software (version 3.5.1) was used to build co-location and co-expressed network between DE-lncRNAs and DEGs.

Statistical Analysis
Significant differences between two groups with the resulting p values of the DEGs and transcripts were adjusted to p-adj (padj < 0.05). Average data of bacterial plate colony count of each sample was transformed to log 10 (cfu/mg), and based on the differences of mean ± SEM of the data, the origin 2017 software was used to map. Chi-square test was used to determine the differences in relative percent survival (RPS) post the Ea challenge between two infected groups of eels.

The RPS, Bacterial Quantity, and Pathological Changes
At 14 days post the challenge (dpc) of Ea, the mortality rates of eels in two groups were 86% and 43% corresponding to 50% RPS in FCIA_Li vs Con_Li (Fig. 1A). X 2 test (chi-square) showed that there was significant difference (p < 0.05) in cumulative mortality between FCIA_Li and Con_Li group from 4 to 14 dpc. Average bacteria load of 10 4.2 cfu and 10 3.7 cfu per mg liver in the eels of Con_Li group, corresponding to average 10 2.1 cfu and 10 2.2 cfu that of FCIA_Li group at 24 and 48 h post challenge (hpc) (Fig. 1B). In the FCIA_Li group, slight edema of hepatocytes and blood coagulation in dilate hepatic vessels were observed at 48 hpc (Fig. 1C). In the liver of uninfected Li group, no colony was showed on the plate at 24-and 48-h sampling, and the hepatocytes were normal as well as red blood cells kept in the blood vessels at 48-h sampling (Fig. 1B, D). However, thrombi formed in the dilated hepatic vessels, celluloid degeneration of vascular wall, and hepatocyte atrophy with eosinophilic cytoplasm in the Con_Li group were observed in the eels at 48 hpc (Fig. 1E).

Quality of RNA-seq
The quality evaluation of sequencing output data was shown in Supplementary Table 1. Compared to the genome DNA of A. anguilla, the average sequencing base number of livers in two infection group was more than 14.5G, but only 9.6 G in the uninfected livers. The average total mapping of 9 samples was more than 81.8 according to the annotation Fig. 1 The cumulative mortality rate of eels challenged by Ea B79 (A), bacteria load in the liver (log 10 cfu/mg) (B), and HE staining results of liver sections of eels in FCIA_Li (C), Li (D), and Con_Li (E) group. A Cumulative mortality rate of eels challenged by Ea B79 in the Li and FCIA groups; *p < 0.05 between the FCIA_Li and the Con_Li group in 14 dpc. B Colony quantity per mg liver tissue of 2 challenged groups and colony growth of representative coating plates of 3 groups. *p < 0.05 and **p < 0.01 between the FCIA_Li and the Con_Li group. C Slight atrophy of hepatocytes and blood coagulation in the hepatic vessels of eels infected by Ea after FCIA injection. D Normal hepatocytes and red blood cells in normal blood vessels of uninfected eels. E Hepatic vein dilation, thrombosis, atrophy and degeneration of hepatocytes, and necrosis of some hepatocytes of eels infected by Ea. Scale bar: 50 μm file of A. anguilla (. gff). In all 9 samples, the base error rate of more than 94.35 fragments was less than 1% (Q20), and the base error rate of more than 87.97% fragments was less than 1‰ (Q30). Mapping rate of 9 samples between reads and genome of A. anguilla were also showed in the Supplementary Table 1.

Heat Map Clusters of Differentially Expressed Genes and Transcripts
From the comparison to the sequencing data of A. anguilla genome, the results of two-way cluster analysis showed quantity of DEGs and DE-mRNAs was about 10 times as the DE-lncRNAs and DE-TUCPs between 3 groups (Fig. 2). The results of DEGs and DE-mRNAs clustering showed the samples in the FCIA_Li group were similar to that of the Li group, but distinct differences were observed between these 6 samples and 3 samples in the Con_Li group ( Fig. 2A, B). Simultaneously, the clustering results of DE-lncRNAs and DE-TUCPs (Fig. 2C, D) showed samples between 3 groups were different contrast to 3 samples in the lncRNAs were similar. Therefore, compared to the results of Gene, mRNA and TUCP clustering, the difference of the lncRNA between inter groups is much greater than that of intra groups.

DE-lncRNAs and the Target DEGs
In total of 9 samples, 13,499 lncRNAs were expressed including 10,176 annotated and 3423 novel lncRNAs (Supplementary Table 2). The Venn map showed 61, 63, and 76 DE-lncRNAs independently expressed in 3 comparisons of Con_Li vs Li, FCIA_Li vs Li, and FCIA_Li vs Con_Li, corresponding to 22, 25, and 28 cross lncR-NAs between the 3 comparisons ( Fig. 3A), and in the volcano map, there were 44, 71, and 34 significantly upregulated corresponding to 67, 58, and 76 downregulated in the 3 comparisons (Fig. 3B). Compared to the Con_Li group, 88 and 96 GO terms were significantly enriched (qvalue < 0.05) in co-expression (Supplementary Table 3) and co-location (Supplementary Table 4) DEGs targeted by 129 DE-lncRNAs in the FCIA_Li group. Top 50 enrichment GO terms of Co-expression showed 10 processes in BP, 4 processes in CC and the protein binding process in MF were distinctively enriched by more than 1000 regulated DEGs (Fig. 3C), and top 52 enrichment GO terms showed protein binding process in MF was distinctively enriched by 346 regulated DEGs in co-location (Fig. 3D). In the comparison of FCIA_Li vs Con_Li, it was showed that 129 DE-lncRNAs targeted 7742 genes in co-expression (Supplementary Table 5), of which 8 KEGG pathways were significantly enriched (Fig. 3E) and the input vs background gene number were higher than 0.376. The co-location of 129 DE-lncRNAs target 942 genes in this comparison (Supplementary Table 6), and the KEGG pathway of cell adhesion molecules (CAMs) was significantly enriched (Fig. 3F).

Validation of qRT-PCR for RNA_seq
In 23 verified DEGs, the fold change of 20 DEGs was kept at the same level, but 2 DEGs (FOS like and serpin peptidase) were higher and the solute carrier family 7 was lower in RNA-seq contrast to the result of qRCR (Fig. 4). Compared to the Con_Li group, a lot of immune-related genes were markedly upregulated in the FCIA_Li group, including complement component, lysozyme C, alpha-2-macroglobulin, H-2 class histocompatibility antigen, interferon regulatory factor, Toll-like receptor, Neutrophil cytosol factor, mitogenactivated protein, and C-C motif chemokine 17 (Fig. 4A, C). The linear fit of 23 genes expression showed correlative coefficient between the two methods was nearly 0.7 (Fig. 4B), and 17 genes were significantly unregulated as well as 4 genes were downregulated in the FCIA_Li group (Fig. 4C).

Discussion
In recent years, studies on the septicemia of cultivated fish caused by Ea have been reported (Armwood et al. 2019;Shao et al. 2015), but the molecular mechanism of the antiinfection in eel is not clear. In a previous study, we found that the Ea mainly invaded the liver of eels in the early stage of infection (Zhai et al. 2021), and the bacteria also caused purulent infection in the liver. Therefore, the liver was selected for RNA-seq research. In this study, the infected liver was found serious pathological changes contrast to eels in the FCIA_Li group showed slight pathological changes due to bacteria load of FCIA_Li group was only 10% as that of the Con_Li group at 24 and 48 hpc (Fig. 1). Based on the RNA-seq of mRNA and lncRNA, the lncRNAs of European eels in the comparison of FCIA_Li vs Con_Li was preliminarily explored, and the anti-infection-related DEGs targeted by DE-lncRNAs were firstly analyzed to reveal the anti-infection effect of the FCA and FIA adjuvant post the challenge of Ea. LncRNA has been found to play an important role in the pathogenesis of many diseases (Chen et al. 2013). From the functional point of view, lncRNA mainly acts on the transcription, post transcription, translation, and apparent level of the target gene respectively through Cis and Trans regulation to regulate the expression of the target gene (Florian and Joshua 2018), so as to find out the different GO terms or KEGG pathways of different treatments through the regulation of lncRNA (Schmittand and Chang 2016). Co-location refers to the fact that lncRNA may regulate the nearby protein coding genes through analyzed the genes within 100-kb upstream and downstream of lncRNA (Bao et al. 2019). In this study, the co-expression and co-location of lncRNA and mRNA were used to predict the target genes of all DE-lncRNAs. The results showed that there were 129 DE-lncRNAs between 2 infection groups of eels, and the target DEGs in co-expression and co-location were as high as 7742 and 942, respectively.
In this study, we analyzed the DEGs, DE-mRNAs, and DE-lncRNAs, and the results showed distinctive difference between Li and FCIA_Li groups in DEGs and DE-mRNAs, especially in DE-lncRNAs (Fig. 2) which was selected to investigate the cause of lower mortality in the FCIA _Li group. So far, there is no report on the molecular mechanism of the protection effect of Freund's adjuvant, and this study will preliminarily reveal the non-specific protective mechanism of the adjuvant.
European eels were inoculated with the FCA and FIA at 0 day and 14 days respectively to boost the protective effect to the Ea infection, and the results showed the RPS was 50% in European eels, much higher than 11% RPS after Ea challenge post 28 d of the solo FIA inoculated (He et al. 2020a, b). Reason for RNA-seq at 48 hpc was that eels start to die from the third day in the FCIA group resulting in unavailable sampling (Fig. 1A). Compared to the Li group, most of the DE-lncRNAs in two infected groups were downregulated (Fig. 3B, 44 and 34 significantly upregulated corresponding to 67 and 76 downregulated). However, volcano map showed 71 upregulated and 58 downregulated lncRNAs in the comparison of FCIA_Li vs Con_Li group. Therefore, the DE-lncRNAs between two infected groups more likely reveal the anti-infection molecular mechanism of the Freund's adjuvant to the Ea challenge (Huang et al. 2015).
The results showed that 88 and 96 GO terms were enriched (padj < 0.05) in the co-expressing and co-location target genes respectively, including cellular response to stimulus, membrane, protein binding, and binding processes et al. (Fig. 3C, D), of which the process of cellular response to stimulus was directly correlate with the immune function of the host (Liu et al. 2017). Eight enriched KEGG pathways were found in the co-expression DEGs targeted by lncRNAs (Fig. 3E), among which the pathways of neuroactive ligandreceptor interaction (Wang et al. 2021), calcium signaling pathway and Wnt signaling pathway (De 2011), MAPK signaling pathway (Wu et al. 2019), CAMs (Krishnan et al. 2019), and cytokine-cytokine receptor interaction (Han et al. 2019) were directly related to non-specific immune clearance function. Cell adhesion molecules (CAMs) are glycol-proteins expressed on the cell surface and play a critical role in a wide array of biologic processes that include hemostasis, the immune response, inflammation, embryogenesis, and development of neuronal tissue (Khodabandehlou et al. 2017). There are four main groups: the integrin family, the immunoglobulin superfamily, selectins, and cadherins (Samanta and Almo 2015). Membrane proteins that mediate immune cell-cell interactions fall into different categories, namely those involved in antigen recognition, costimulation, and cellular adhesion (Gennarini and Furley 2017). There is only the CAMs KEGG pathway was enriched in co-location (Fig. 3F) in this study, and 15 DEGs in CAMs targeted by 6 DE-lncRNAs mainly encodes myelinassociated glycoprotein (Bornhöfft et al. 2018), B-cell receptor CD22 , class I histocompatibility antigen (He et al. 2020a, b), and hemicentin proteins (Schultz et al. 2005). These proteins played important role in the process of the host anti-Ea infection.
In this study, 23 immune-related DEGs were randomly selected for qRT-PCR verification. Although 20 DEGs were basically kept at the same level, there were some differences between the results of qPCR verification and RNAseq (Fig. 4) and the fit linear result showed their correlative coefficient is relatively low, which probably due to the multiple transcripts of some gene in eukaryotic cells. Primers used in qPCR is actually designed based on the level of transcripts, but RNA-seq can be used for differential analysis at both gene and transcripts level. For the same gene, qPCR is the difference in the expression level of a single transcript between samples, while RNA-seq is the comprehensive result of the expression level of multiple transcripts due to alternative splicing of a gene. Compared to the Con_Li group, 9 genes of 23 immune-related genes increased more than 10 times in the FCIA_Li group (Fig. 4C), indicating these genes involved in non-specific immunity in eels (Zhao et al. 2020;He et al. 2020a, b;Xiao et al. 2022). The interaction between DE-lncRNAs of 2 infection groups and the target DEGs through co-expression and colocation showed 129 DE-lncRNAs-targeted 8684 DEGs. To compare the possible anti-bacterial function of the target DEGs, we selected 16 DE-lncRNAs changed more than 20 times between 2 infected groups (Figs. 5 and 6; Table 2). The result showed most of DEGs targeted by upregulated DE-lncRNAs in the FCIA_Li group were directly or indirectly involved in the immune function of the host, but that of DEGs targeted by upregulated DE-lncRNAs in the Con_Li group were closely related to immunosuppression or apoptosis. Some of the upregulated DEGs were also verified by the qRT-PCR method. Therefore, the Freund's adjuvant enhanced the non-specific immune function of European eel and prevented the apoptosis of hepatocytes after bacterial infection, but immune function of eels in the control group was suppressed after the Ea infection and hepatocytes were apoptotic due to pro-inflammatory factors, which eventually led to the increased mortality of infected eels.

Conclusion
Compared to the eels inoculated by FCA and FIA, thrombi in the hepatic vessels, necrosis, and atrophy of hepatocytes were observed, and a large number of bacteria were isolated from the control eels corresponding to 50% RPS in FCIA_ Li vs Con_Li. The strand-Specific RNA-seq showed obvious difference appeared in the clustering of DE-lncRNAs contrast to that of DEGs and DE-mRNAs. In total, 13,499 lncRNAs were expressed in 9 samples with 10,176 annotated and 3423 novel lncRNAs, of which 129 DE-lncRNAs were highlighted between two infected groups of FCIA_Li and Con_Li. One hundred and eighty-four enriched GO terms indicated DEGs targeted by 129 DE-lncRNAs in co-expression and co-location were mainly involved in single-organism cellular process, biological regulation, regulation of cellular process and response to stimulus in BP, component of membrane in CC, and protein binding in MF. The target DEGs were enriched in Eight KEGG pathways in co-expression and enriched in cell adhesion molecules pathway in co-location. The interaction between 129 DE-lncRNAs and 8684 target genes in the comparison of FCIA_Li vs Con_Li showed 118 DE-lncRNAs target 1161 DEGs and formed a total of 8474 and 333 related links in co-expression and co-location network, of which 16 crucial DE-lncRNAs linked with 860 DEGs in co-expression and 79 DEGs in co-location were further analyzed. The DE-lncRNAs revealed in our study will accelerate the understanding of the underlying Edwardsiella infection and fish anti-infection processes.

Ethics Approval
The fish used in this study were cultivated under protocols approved by the Jimei University Animal Care Committee , and all experiments complied with the current laws of China and international organizations for the use of experimental animals.

Competing Interests
The authors declare no competing interests.