Multi-Organ Transcriptome Response of Lumpfish (Cyclopterus lumpus) to Aeromonas salmonicida Subspecies salmonicida Systemic Infection

Lumpfish is utilized as a cleaner fish to biocontrol sealice infestations in Atlantic salmon farms. Aeromonas salmonicida, a Gram-negative facultative intracellular pathogen, is the causative agent of furunculosis in several fish species, including lumpfish. In this study, lumpfish were intraperitoneally injected with different doses of A. salmonicida to calculate the LD50. Samples of blood, head-kidney, spleen, and liver were collected at different time points to determine the infection kinetics. We determined that A. salmonicida LD50 is 102 CFU per dose. We found that the lumpfish head-kidney is the primary target organ of A. salmonicida. Triplicate biological samples were collected from head-kidney, spleen, and liver pre-infection and at 3- and 10-days post-infection for RNA-sequencing. The reference genome-guided transcriptome assembly resulted in 6246 differentially expressed genes. The de novo assembly resulted in 403,204 transcripts, which added 1307 novel genes not identified by the reference genome-guided transcriptome. Differential gene expression and gene ontology enrichment analyses suggested that A. salmonicida induces lethal infection in lumpfish by uncontrolled and detrimental blood coagulation, complement activation, inflammation, DNA damage, suppression of the adaptive immune system, and prevention of cytoskeleton formation.

Bacterial diseases are a health concern for lumpfish delousing performance and survival rates in sea-cages [3,4]. Aeromonas salmonicida, a facultative intracellular pathogen, endemic worldwide in fresh and marine water, and the etiologic agent of furunculosis [7][8][9][10], incubated at 15 • C for 4 days to determine the number of A. salmonicida CFU per g of tissue or 1 mL of blood, respectively. The Tukey's multiple comparisons test followed one-way ANOVA (p ≤ 0.05 was considered significant). Statistical analyses and data visualization were performed using GraphPad Prism 7.0 (GraphPad Software, La Jolla, CA, USA).

Histopathology
Sections of lumpfish head-kidney, spleen, and liver were collected from the noninfected control and infected fish at 3 and 10 dpi. Tissues were submerged in 10% phosphate saline buffered (PBS) formalin for three consecutive days at room temperature. The fixed tissues were removed from formalin and stored in PBS at 4 • C until processing for paraffinembedded tissue block according to established procedures [41]. Tissues were sectioned, and 5 µm sections were stained with hematoxylin and eosin (Leica Biosystems, Ontario, Canada) and visualized under a light microscope (Olympus CX21, New York, NY, USA).

RNA Purification
For RNA-sequencing (RNA-Seq) and real-time quantitative polymerase chain reaction (qPCR) analyses, the tissue (~100 mg) of head-kidney, spleen, and liver were sampled from three individual lumpfish (non-infected fish), and at 3, and 10 dpi with 10 4 CFU dose −1 , similar to other studies [42][43][44][45][46][47]. Lumpfish tissues (n = 27) were preserved in 500 µL of RNAlater according to the manufacturer's instructions (Invitrogen, Carlsbad, CA, USA) until further processing. RNA was extracted from fish tissues using the mirVana RNA isolation kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions. RNA samples were treated with 2 U mL −1 of TURBO DNase (TURBO DNA-free™ Kit, Invitrogen, Carlsbad, CA, USA) following the manufacturer's instructions for the complete digestion of DNA and removal of remaining DNase and divalent cations, such as magnesium and calcium. Purified RNA samples were quantified for concentration and evaluated for purity using a spectrophotometer (Genova-nano, Jenway, Stone, Staffordshire, England), and evaluated for integrity using 1% agarose gel electrophoresis (Table S1).

Library Preparation and RNA-Seq
For each group, three biological replicates were analyzed for a total of 27 samples ( Figure S1). Library construction and sequencing services were performed at Genome Quebec, Quebec, Canada. Briefly, high-quality RNA was evaluated using a Bioanalyzer 2100 (Agilent, Santa Clara, CA, USA), and only samples with a minimum RIN of 5 proceeded to the library construction ( Figure S2, Table S1). Libraries were generated from 250 ng of total RNA using the NEBNext ® Multiplex Oligos for Illumina ® (Dual Index Primers Set 1; Adapter 1: 3'-AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC-5'; Adapter 2: 3'-AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT-5') and sequenced in a NovaSeq 6000 Sequencer (Illumina) using a NovaSeq 6000 S4 100 bp PE flow cell. The raw data were deposited in the NCBI Sequence Read Archive (SRA) (accession number PRJNA596743).

Reference Transcriptome Assembly and Downstream Analysis
After RNA sequencing, paired-end raw reads were mate-paired and filtered to remove low-quality reads using the CLC Genomics Workbench v20.0 (CLCGWB; Qiagen, Hilden, Germany) with default parameters (mate-paired read information, minimum distance = 1; maximum distance = 1000). Adapter trimming was done by the CLCGWB using the trim reads tool with default parameters (quality trimming, trim using quality scores, limit: 0.05, and trim ambiguous nucleotides, the maximum number of ambiguities = 2). The quality of the reads was checked using the FastQC [48], and a multiqc [49] report was generated before and after trimming. CLCGWB then mapped trimmed reads against the lumpfish genome (NCBI accession number PRJNA625538) using the RNA-Seq analysis tool. Reads mapping and transcript counts were conducted using the following settings: mismatch cost = 2, insertion and deletion costs = 3, minimum length fraction and minimum similarity fraction = 0.8, the maximum number of hits for a read = 10, and strand-specific = both. Gene Microorganisms 2022, 10, 2113 5 of 32 expression quantification and normalization of the mapped reads were performed by an alignment-dependent expectation-maximization (EM) algorithm [50] based on the RSEM and eXpress methods [51]. The TPM values were then computed from the counts assigned to each transcript after normalization by the trimmed mean of M-values (TMM) [52]. A global correlation analysis was performed using log 2 -transformed TPM values (x + 1) of each gene, and transcripts and correlations were estimated by the Pearson method. Abundance data were subsequently subjected to differential expression analyses using the CLCGWB and the differential expression tool based on a negative binomial general linear model (GLM) [53]. A standard selection of biologically significant differentially expressed genes (DEGs) and differentially expressed transcripts (DETs) were performed with cut-off values of log2 fold-change (FC) ≥ 1 and a false discovery rate (FDR) p ≤ 0.05.

De Novo Transcriptome Assembly, Contig Abundance, and Functional Annotation
Adaptor sequences and reads ≤50 base pairs (bp) were trimmed using the Trimmomatic v0.38 [54]. Reads were assembled into a de novo transcriptome using the Trinity software (v2.8.4) with default parameters [55]. Assembled contigs shorter than 200 bp were excluded from the analysis. We examined the read representation of the assembly by aligning the processed reads on de novo assembled transcriptome using the samtools v1.9 [56] and bowtie2 v2.3.5 [57]. We also inspected the read supports for assembled transcripts using IGV v2.7.2 [58]. Furthermore, we also examined how successfully assembled proteincoding transcripts were reconstructed to full-or near full-length using the BLAST+ and Swissprot/TrEMBL. The de novo transcriptome assembly quality and completeness were evaluated using BUSCO version 3 [59] against a predefined set of 4584 Actinopterygian single copy orthologs from the orthoDB v9 database. RSEM (v1.3.1) [60] was used to quantify TPMs in the trinity package. Trinotate v3.1.1, a functional annotation pipeline, was used to generate an annotation report for the potential biological function of the assembled contigs [61]. Trinotate uses the TransDecoder v5.5 [55] to identify protein-coding regions in each assembled transcript. BLAST+ variants (blastn, blastx, and blastp) against sequence databases downloaded locally (last accessed, January 2021), including RefSeq-rna, RefSeq-protein, nt, nr and SwissProt, were used to annotate the de novo assembled transcriptome. The transcriptomes were further annotated for remote homologs and protein domains using HMMER v3.2.1 and Pfam v3.2.1 [62,63]. The SignalP v5.0 [64] and tmhmm v2.0c [65] software tools were used to predict signal peptides and transmembrane domains, respectively.
A differential expression (DE) analysis was performed using DESeq2 [66] and edgeR [53]. Within each pairwise comparison, only transcripts with an FDR adjusted p-value ≤ 0.05 were considered significantly differentially expressed. The DE analysis was conducted at the isoform level to identify the DETs in each organ at different timepoint.

Gene Filtration and Gene Ontology (GO) Enrichment Analysis
DEGs identified by reference transcriptomic assembly analysis were filtered (log 2 FC ≥ 1, p-value ≤ 0.05) to identify the enriched Gene Ontology terms. On the other hand, DETs identified by de novo reference assembly analysis and edgeR analysis were utilized for the nucleotide blast against the lumpfish genome in NCBI to extract the corresponding gene symbols. Those genes were added to the Gene Ontology term enrichment analyses. To obtain an overall view, GO enrichment analysis of all DEGs of head-kidney, spleen, and liver at 3 dpi and 10 dpi were conducted by the ClueGO App [67] (2.5.8) in Cytoscape 3.9 [68] using ClueGO source files for lumpfish (Supplementary Data 1). ClueGO source files were created using a GO OBO file downloaded on 24 March 2022. Fisher's exact test was conducted to study the enrichment of GO terms with a p-value cut off ≤0.05 for 3 dpi and a p-value cut off ≤0.00001 for 10 dpi. A differential p-value cut-off and the GO term fusion strategy were employed to reduce the redundancy of the GO terms and simplify the network. To obtain a global view, GO enrichment analysis of DEGs (n = 600 for 10 dpi; 300 most significant DEGs (lowest FDR p-value) from each assembly analysis) of Microorganisms 2022, 10, 2113 6 of 32 head-kidney, spleen, and liver at 3 dpi and 10 dpi were conducted by the ClueGO. Fisher's exact test was conducted to study the enrichment of GO terms with a p-value cut off ≤0.05. The GO term fusion strategy was employed. However, to explore the pathogenesis in-depth, GO enrichment analysis associated with biological processes of upregulated and downregulated DEGs in the head-kidney 3 dpi and 10 dpi, spleen 3 dpi and 10 dpi, and liver 3 dpi and 10 dpi was conducted by setting the network specificity at medium in the ClueGO. Fisher's exact test was conducted to study the enrichment of GO terms with a p-value cut off ≤0.05 (3 dpi) and 0.001 (10 dpi). The GO term fusion strategy was applied.

Real-Time Quantitative Polymerase Chain Reaction (qPCR) Analyses
To verify the RNA-Seq analyses, expression levels of 14 genes were measured in the same 27 RNA samples that were subjected to RNA-Seq analyses. These genes were selected based on their TPM values as they were expressed in individual samples from at least one group (i.e., uninfected control, 3 dpi or 10 dpi) in all three tissues (i.e., head-kidney, liver, and spleen). In all cases, first-strand cDNA templates were synthesized, and qPCR amplifications were performed as described previously [20,23].
The sequences, amplicon sizes, and efficiencies for all primer pairs used in the qPCR analyses are presented in Table S2. Primer pairs for il8b, tlr5a, saa5/app, and the endogenous control transcripts were designed, and quality control (QC) was tested previously [20]. Primers new to this study were designed following the same parameters. All primer pairs used herein were (re)-subjected to QC testing [20] using a cDNA pool generated from the spleen control group and one from the spleen 10 dpi group. The calculated efficiencies are an average of the two values.
In the experimental qPCR analyses, expression levels of the genes were normalized to expression levels of two endogenous control transcripts. The fluorescence threshold cycle (C T ) values of all 27 samples in the study were measured (in triplicate) for each of these transcripts using cDNA of 4 ng of input total RNA, and then analyzed using geNorm [69]. Based on this analysis, eukaryotic translation initiation factor 3 subunit D (etf3d) (geNorm M = 0.31) and elongation factor 1-alpha (ef1a) (geNorm M = 0.34) were selected as the two endogenous controls.
Primer QC and endogenous control selection were followed by the experimental qPCR analyses. cDNA representing 4 ng of input RNA was used as a template in the PCR reactions. On each plate, for every sample, the selected genes and endogenous controls were tested in triplicate, and a non-template control was included. The relative quantity (RQ) of each transcript was determined using the QuantStudio Real-Time PCR Software (version 1.3) (Applied Biosystems), where Ct values were normalized with both etf3d and ef1a with amplification efficiencies incorporation. For each target of interest (TOI), the sample with the lowest normalized expression (mRNA) level was set as the calibrator sample (i.e., assigned an RQ value = 1).
To compare the TPM to the RQ values, the log 2 normalized values were utilized. Statistical regression analyses and data visualization were performed using the GraphPad Prism 7.0 (GraphPad Software, La Jolla, CA, USA).

LD 50 Determination and A. salmonicida Infection Kinetics in Lumpfish
Five groups of 60 fish (duplicated tanks, total n = 600) were injected with five different doses of A. salmonicida J223 (10 1 , 10 2 , 10 3 , 10 4 , and 10 5 CFU dose −1 ) to determine the LD 50 and infection kinetics ( Figure 1A). After 3 days post-infection (dpi), lack of appetite, erratic swimming, and internal hemorrhagic septicemia were observed. Fish started to die at 7-10 dpi, and there was no survival in the fish infected with the 10 3 , 10 4 , and 10 5 CFU dose −1 . 32% of lumpfish survived after the infection with 10 2 CFU dose −1 ( Figure 1A). In contrast, 93% of fish survived after the infection with the lowest dose tested (10 1 CFU dose −1 ) ( Figure 1A). The survivors and the non-infected control fish showed no symptoms of disease or mortality. The LD 50 for A. salmonicida J223 in lumpfish was 187 CFU dose −1 according to the Reed and Muench method [39] and 273 CFU dose −1 according to the Karber method [40]. According to the log-rank (Mantel-Cox) test and log-rank test for trend, the survival curves and the trend were significantly different (p < 0.0001). and infection kinetics ( Figure 1A). After 3 days post-infection (dpi), lack of appetite, erratic swimming, and internal hemorrhagic septicemia were observed. Fish started to die at 7-10 dpi, and there was no survival in the fish infected with the 10 3 , 10 4 , and 10 5 CFU dose −1 . 32% of lumpfish survived after the infection with 10 2 CFU dose −1 ( Figure 1A). In contrast, 93% of fish survived after the infection with the lowest dose tested (10 1 CFU dose −1 ) (Figure 1A). The survivors and the non-infected control fish showed no symptoms of disease or mortality. The LD50 for A. salmonicida J223 in lumpfish was 187 CFU dose −1 according to the Reed and Muench method [24] and 273 CFU dose −1 according to the Karber method [25]. According to the log-rank (Mantel-Cox) test and log-rank test for trend, the survival curves and the trend were significantly different (p < 0.0001).  Lumpfish were ip injected with five different doses of A. salmonicida ranging from 1.1 × 10 1 to 1.1 × 10 5 CFU dose −1 . The survival percentage of lumpfish infected with these different doses were 93%, 32%, 0%, 0%, and 0% at 30 dpi, respectively; (B-E). A. salmonicida colonization in the lumpfish lymphoid organs (head-kidney, spleen, and liver) and blood at 3, 7, 10, 21, 33 dpi; (B). A. salmonicida colonization in head-kidney; (C). A. salmonicida colonization in spleen; (D). A. salmonicida colonization in liver; (E). A. salmonicida colonization in blood; *: indicates a significant statistical difference (p < 0.05); (F). Histopathology of lumpfish tissues stained with hematoxylin and eosin. Lumpfish tissues were collected from non-infected control fish and infected fish (1.1 × 10 4 CFU dose −1 ) at 10 dpi, visualized under the light microscope (× 400). Blue and red arrows indicate inflammatory cells and necrosis, respectively.
A. salmonicida colonization was determined at 0, 3, 7, 10, 14, 21, and 33 dpi in different tissues ( Figure 1B-E). A. salmonicida was detected in the head-kidney at 3 dpi in fish infected with 10 3 -10 5 CFU dose −1 , but not when infected with 10 1 -10 2 CFU dose −1 . We did not detect bacteria in the liver, spleen, and blood at 3 dpi in lumpfish injected with the lowest doses (10 1 -10 4 CFU dose −1 ) ( Figure 1C-E). In contrast, a few fish infected with the 10 5 CFU dose −1 showed bacterial colonization in the spleen and liver at 3 dpi ( Figure 1C,D). A. salmonicida was first detected in blood at 7 dpi ( Figure 1E). At 7 dpi, head-kidney, spleen, and liver showed bacterial colonization in all doses tested except in the lowest doses evaluated (10 1 and 10 2 CFU dose −1 ) ( Figure 1B-D). At 10 dpi, bacterial colonization was detected in all the tissues, except for the fish infected with 10 1 CFU dose −1 ( Figure 1B-E). Significant differences between bacterial loads were observed only in the spleen samples of the fish infected with 10 5 CFU dose −1 compared to the other groups (p < 0.003). At 14 dpi, 3 fish infected with the 10 3 CFU dose −1 showed bacterial colonization in all tissues sampled ( Figure 1B-E). After 15 dpi, no bacteria were detected in the remaining survivors in the 10 2 CFU dose −1 infected group. However, A. salmonicida was detected in the head-kidney, spleen, and liver until 33 dpi in fish infected with 10 1 CFU dose −1 . These results suggest that A. salmonicida targets the head-kidney and disseminates to other organs, causing a systemic infection.
Histopathological analysis indicated that A. salmonicida caused inflammation and tissue necrosis in head-kidney, spleen, and liver ( Figure 1F). Previously, we observed intracellular A. salmonicida in lumpfish tissues at 3 and 10 dpi [70], and similar observations were determined in the current study.

Raw Sequencing Data and Quality Statistics
RNA samples were collected from the head-kidney, spleen, and liver of three noninfected lumpfish and three infected lumpfish (10 4 CFU dose −1 ) at 3 and 10 dpi ( Figure S1). RNA quality is shown in Supplementary Files (Table S1 and Figure S2A-C). RNA sequencing generated 1.08 billion Illumina NovaSeq reads ranging from 67-95 million raw reads per sample (Table S3), with a length of 101 bp. After trimming, reads were subjected to reference genome-based transcriptome assembly analysis and de novo transcriptome assembly analysis ( Figure S1).

Global Profile of Differentially Expressed Genes and Transcripts Identified Using the Lumpfish Reference Genome
To study the lumpfish response to A. salmonicida infection, we profiled the global gene expression of the head-kidney, spleen, and liver at 3 and 10 dpi, compared to non-infected fish using RNA-Seq. A global gene expression correlation analysis showed a high degree of correlation (R 2 = 0.89 to 0.98; p < 0.0001) between different experimental conditions ( Figure S3). Principal component analysis (PCA) and heatmap results reveal a clear tissue and time point clusterization (Figures 2 and 3). Among all the organs studied, the spleen has the clearest clusterization ( Figure 3). The log2 fold-change (FC) ≥ 1 and false discovery rate (FDR) p-value of ≤0.05 were set as the cut-off criteria for sorting out significant differentially expressed genes. We found 102 differentially expressed genes (DEGs) in the head-kidney at 3 dpi. These DEGs included 94 upregulated and 8 down-regulated genes (Table 1, Figure 4A). Also, 1922 DEGs were identified in the head-kidney at 10 dpi. These DEGs included 530 upregulated and 1392 downregulated genes (Table 1, Figure 4A). In the spleen, 637 DEGs were identified at 3 dpi, including 253 upregulated and 384 downregulated genes (Table 1, Figure 4B). In  (C). PCA of lumpfish spleen infected with A. salmonicida, green dot represents three control samples, red dot represents three 3 dpi samples, blue dot represents three 10 dpi samples; (D). Heat map of lumpfish spleen infected with A. salmonicida; (E). PCA of lumpfish liver infected with A. salmonicida, green dot represents three control samples, red dot represents three 3 dpi samples, blue dot represents three 10 dpi samples; (F). Heat map of lumpfish liver and liver infected with A. salmonicida.

Figure 3. Clusterization of gene expression profile in different lumpfish organs infected A. salmonicida. (A). Principal component analysis (PCA) of lumpfish head-kidney infected with
A. salmonicida, green dot represents three control samples, red dot represents three 3 dpi samples, blue dot represents three 10 dpi samples; (B). Heatmap of lumpfish head-kidney infected with A. salmonicida; (C). PCA of lumpfish spleen infected with A. salmonicida, green dot represents three control samples, red dot represents three 3 dpi samples, blue dot represents three 10 dpi samples; (D). Heat map of lumpfish spleen infected with A. salmonicida; (E). PCA of lumpfish liver infected with A. salmonicida, green dot represents three control samples, red dot represents three 3 dpi samples, blue dot represents three 10 dpi samples; (F). Heat map of lumpfish liver and liver infected with A. salmonicida.
The log 2 fold-change (FC) ≥ 1 and false discovery rate (FDR) p-value of ≤0.05 were set as the cut-off criteria for sorting out significant differentially expressed genes. We found 102 differentially expressed genes (DEGs) in the head-kidney at 3 dpi. These DEGs included 94 upregulated and 8 down-regulated genes (Table 1, Figure 4A). Also, 1922 DEGs were identified in the head-kidney at 10 dpi. These DEGs included 530 upregulated and 1392 downregulated genes (Table 1, Figure 4A). In the spleen, 637 DEGs were identified at 3 dpi, including 253 upregulated and 384 downregulated genes (Table 1, Figure 4B). In the spleen, 3133 DEGs were identified at 10 dpi, including1368 upregulated and 1765 downregulated genes (Table 1, Figure 4B). In the liver, 58 DEGs were identified at 3 dpi. These DEGs included 44 upregulated and 14 downregulated genes (Table 1, Figure 4C). Also, 2766 DEGs were identified in the liver at 10 dpi, including 1360 upregulated and 1406 downregulated genes (Table 1, Figure 4C). Gene identifier, description/annotation, fold-change, and FDR (p-value) are listed in File S1. A comparison between the head-kidney, spleen, and liver, including 3-and 10-dpi, showed that 309 DEGs were common to all organs, while the head-kidney and the spleen shared 373 DEGs, the head-kidney and the liver shared 196 DEGs, and the spleen and the liver shared 738 DEGs ( Figure 4D).  Similarly, the log2FC ≥ 1 and FDR p-value of ≤0.05 were set as the cut-off criteria for sorting out significant differentially expressed transcripts (DETs). We identified 133 DETs in the head-kidney at 3 dpi, including 89 upregulated and 44 down-regulated transcripts ( Table 2). Also, 699 DETs were identified in the head-kidney at 10 dpi, including 451 up- Similarly, the log 2 FC ≥ 1 and FDR p-value of ≤0.05 were set as the cut-off criteria for sorting out significant differentially expressed transcripts (DETs). We identified 133 DETs in the head-kidney at 3 dpi, including 89 upregulated and 44 down-regulated transcripts (Table 2). Also, 699 DETs were identified in the head-kidney at 10 dpi, including 451 upregulated and 248 downregulated transcripts ( Table 2). In the spleen, 614 DETs were identified at 3 dpi, including 237 upregulated and 377 downregulated transcripts ( Table 2). In the head-kidney, 2152 DETs were identified at 10 dpi, including 1173 upregulated and 979 down-regulated transcripts ( Table 2). In the liver, 33 DETs were identified at 3 dpi, including 27 upregulated and six downregulated genes (Table 2). Also, 1697 DETs were identified in the liver at 10 dpi, including 887 upregulated and 810 downregulated transcripts ( Table 2). Gene identifier, fold-change, and FDR (p-value) are listed in File S2.

Global Profile of Differentially Expressed Transcripts Identified by De Novo Transcriptome
To identify potential novel genes and transcripts, a de novo transcriptome analysis was conducted. Quality filtering and trimming were performed by trimmomatic, and approximately 4.28% of the raw reads were removed (Table S4). The remaining highquality reads (originating from the three different lymphoid tissues) were used to build a de novo transcriptome assembly using Trinity v2.8.4 assembler.
The de novo assembly resulted in 403,204 transcripts with an average read length of 497 bp, representing 270,150 genes identified by Trinity (Table S4). The total length of all assembled transcripts is 522,614,427 bp with an N50 length of 3235 bp and GC content of 45.99%. We found that more than 98% of the reads were successfully aligned consistently for each sample (Table S5). Coding transcripts assessment was performed using the blastx search program in the database NCBI, RefSeq RNA, and SwissProt [71,72].
We further evaluated the completeness of the transcriptome assembly using BUSCO. Busco pipeline for gene set completeness was assessed for eukaryotes (n = 303), vertebrates (n = 2586), and actinopterygian (n = 4584). The analysis reported that the majority of the actinopterygian core genes had been successfully recovered from the lumpfish de novo assembly. Specifically, of the 4584 single-copy orthologs searched,~88% were completely recovered, and~4% were partially recovered. Only~8% of single-copy orthologs were classified as missing in the assembly. This data indicates a complete, consistent, high-quality lumpfish transcriptome assembly ( Figure S4).
DETs identified from the three lymphoid tissues at different time points are summarized in File S3. The lists of DETs identified by DESeq2 were generally higher and were almost accommodated within the edgeR DETs lists. We used the more conservative edgeR-generated DETs for further analysis ( Table 3). The log 2 fold-change (FC) ≥ 1 and FDR p-value of ≤0.05 were set as the cut-off criteria for sorting out significant differentially expressed transcripts. We found 286 DETs in the head-kidney 3 dpi, including 138 upregulated and 148 downregulated transcripts (Table 3, Figure S5A). Also, 477 DETs were identified in the headkidney at 10 dpi, including 204 upregulated and 273 down-regulated transcripts ( Table 3, Figure S5A). In the spleen, 501 DETs were identified at 3 dpi, including 214 upregulated and 287 down-regulated transcripts (Table 3, Figure S5B). Also in the spleen, 2415 DETs were identified at 10 dpi, including 1005 upregulated and 1410 downregulated transcripts (Table 3, Figure S5B). In the liver, 133 DETs were identified at 3 dpi, including 56 upregulated and 77 downregulated transcripts (Table 3, Figure S5C). Also in the liver, 2093 DETs were identified at 10 dpi, including 1053 upregulated and 1040 downregulated transcripts (Table 3, Figure S5C). A comparison between all-time points of the head-kidney, spleen, and liver showed 56 DETs in common, while the head-kidney and spleen shared 53 DETs, the head-kidney and liver shared 26 DETs, and the spleen and liver shared 274 DETs ( Figure S5D).
The hierarchal cluster of DETs expressed in abundance (log 2 FC ≥ ±5 and FDR ≤ 0.05) visualized as in the heatmap supports the tissue and time point specific clustering ( Figure S6A). The heatmap also reveals that samples from uninfected lumpfish and infected animals clustered mostly within each tissue sub-cluster ( Figure S6A). Also, we observed that the transcriptomic response was clearly separated based on infection time points in the spleen ( Figure S6A). On the other hand, transcript responses in head-kidney and liver samples from the pre-infected fish and infected fish at 3 dpi were not highly differentiated, indicating an early process of infection in these tissues ( Figure S6A). We also assessed and visualized inter and intragroup variability using Pearson's correlation plots of correlation values between samples that agree with the hierarchal clustering analysis ( Figure S6B).
Furthermore, a blastn analysis of all DETs identified by de novo assembly was conducted against the lumpfish genome using Blast+ 2.12.0 to retrieve lumpfish gene symbols corresponding to those transcripts (File S4). The analysis identified 1954 genes that were common to the DEGs identified by the reference genome-guided transcriptome analysis. In total, 1307 unique genes were identified, which included 477 genes in the head-kidney, 825 genes in the spleen, and 679 genes in the liver ( Figure 5 and File S5). These unique genes were added to the DEGs list generated by the reference genome-guided transcriptome for GO enrichment analysis. ducted against the lumpfish genome using Blast+ 2.12.0 to retrieve lumpfish gene symbols corresponding to those transcripts (File S4). The analysis identified 1954 genes that were common to the DEGs identified by the reference genome-guided transcriptome analysis. In total, 1307 unique genes were identified, which included 477 genes in the head-kidney, 825 genes in the spleen, and 679 genes in the liver ( Figure 5 and File S5). These unique genes were added to the DEGs list generated by the reference genome-guided transcriptome for GO enrichment analysis. .

Gene Ontology Enrichment Analysis
Overall, the GO enrichment analysis using a combination of all DEGs at 3 dpi and 10 dpi identified multiple enriched GO terms related to biological process (BP), molecular function (MF), and cellular component (CC) (Figure S7 and S8, and File S6). This result suggests that nucleic acid metabolism and immune responses are mostly affected at the early point of infection. On the other hand, a lethal A. salmonicida infection could modulate lumpfish adaptive immune responses and metabolic processes.
The GO enrichment analysis using all DEGs in the head-kidney at 3 dpi identified GO terms associated with BP (e.g., response to stimulus), MF (e.g., hydrolase activity), and CC (e.g., intracellular anatomical structure) ( Figure 6A, File S6). The upregulated DEGs in the head-kidney at 3 dpi were associated with acute phase response, inflammatory response, complement activation, negative regulation of immune effector process, fibrin clot formation and others (File S7). However, no GO terms were enriched by the downregulated DEGs of the head-kidney at 3 dpi. All the DEGs of the spleen at 3 dpi

Gene Ontology Enrichment Analysis
Overall, the GO enrichment analysis using a combination of all DEGs at 3 dpi and 10 dpi identified multiple enriched GO terms related to biological process (BP), molecular function (MF), and cellular component (CC) (Figures S7 and S8, and File S6). This result suggests that nucleic acid metabolism and immune responses are mostly affected at the early point of infection. On the other hand, a lethal A. salmonicida infection could modulate lumpfish adaptive immune responses and metabolic processes.
The GO enrichment analysis using all DEGs in the head-kidney at 3 dpi identified GO terms associated with BP (e.g., response to stimulus), MF (e.g., hydrolase activity), and CC (e.g., intracellular anatomical structure) ( Figure 6A, File S6). The upregulated DEGs in the head-kidney at 3 dpi were associated with acute phase response, inflammatory response, complement activation, negative regulation of immune effector process, fibrin clot formation and others (File S7). However, no GO terms were enriched by the downregulated DEGs of the head-kidney at 3 dpi. All the DEGs of the spleen at 3 dpi showed their association with several enriched GO terms related to BP (e.g., cell adhesion, defense response, nucleic acid metabolic process), MF (DNA binding), and CC (e.g., extracellular region, external encapsulating structure) ( Figure 6C, File S6). Upregulated DEGs were mostly associated with acute phase response, complement component activation, humoral immune response, inflammatory responses, and many others (File S7). Downregulated DEGs of the spleen at 3 dpi were associated with ribosome assembly, cytoplasmic translation, and oxygen carrier (File S7). Furthermore, the DEGs of the liver at 3 dpi were only associated with response to stress ( Figure 6E and File S6). Upregulated DEGs showed three enriched GO terms, such as acute phase response, chemoattractant activity, and cellular response to interleukin-1 (File S7), and downregulated DEGs of the liver at 3 dpi were not associated with any GO terms.  Furthermore, the 600 most significant DEGs (lowest FDR p-value) in the head-kidney at 10 dpi were associated with eight enriched GO terms related to BP (e.g., defense response, nucleic acid metabolic process), 1 GO term related to MF (nucleic acid binding), and 10 GO terms related to CC (e.g., organelle, nucleoplasm, extracellular region) ( Figure 6B and File S6). The 600 most significant DEGs (lowest FDR p-value) of the spleen at 10 dpi were associated with 16 enriched GO terms related to BP (e.g., immune response, inflammatory response, defense response), four GO terms associated with MF (signaling receptor binding, organic cyclic compound binding, heterocyclic compound binding, and nucleic acid binding), and nine GO terms associated with CC (e.g., organelle, extracellular region) ( Figure 6D, File S6). The 600 most significant DEGs of the liver at 10 dpi were associated with 13 enriched GO terms related to BP (e.g., defense response, small molecule metabolic process, lipid metabolic process, and nucleic acid metabolic process), three GO terms related to MF (catalytic activity, nucleic acid binding, and oxidoreductase activity), and five GO terms related to CC (e.g., organelle, cell junction, extracellular region) ( Figure 6F and File S6).
Furthermore, our GO enrichment analysis indicates that upregulated DEGs of the head-kidney at 10 dpi were associated with acute phase response, complement activation, inflammation, regulation of the apoptosis process, and negative regulation of the immune effector process (File S7). Upregulated DEGs in the spleen at 10 dpi were associated with complement activation, regulation of the apoptotic process, acute-phase response, blood coagulation, and inflammatory response (File S7). Upregulated DEGs of the liver at 10 dpi were associated with acute-phase response and inflammatory response (File S7).
Downregulated DEGs of the head-kidney at 10 dpi were associated with metabolic processes, ion transport, and microtubule bundle formation (File S7). Downregulated DEGs in the spleen at 10 dpi were associated with cytoskeleton organization, nucleic acid metabolic process, ribosome biosynthesis, and translation (File S7). Downregulated DEGs of the liver at 10 dpi were associated with metabolic processes, such as lipid, organic acid, amino acid, DNA, and RNA metabolic process, ion transport, DNA repair, double-strand break repair, and cell cycle (File S7).

Analysis of the Most Significant DEGs
We identified the 300 most significant DEGs based on the lowest FDR p-value in lumpfish head-kidney, spleen and liver (File S9). Our results indicate that the most significantly overexpressed genes in the head-kidney, spleen, and liver were il1b, il8, il10, il6, hamp, haptoglobin (hp), ptx3, collagenase (mmp13b), c7b, and app ( Figure 7). In addition to this, the top significant upregulated genes in the head-kidney at 3 dpi were fibrinogen beta and gamma chain (fbb and fbg) and complement factor B (cfb) (Figure 8). The top significant upregulated gene in the spleen at 3 dpi was tubulin alpha-1A chain (tuba1a), and at 10 dpi were tlr5, coagulation factor IIIa (f3), and socs3a ( Figure 9). The top significant upregulated genes in the liver at 10 dpi were adenosine receptor A3 (adora3) and carcinoembryonic antigen-related cell adhesion molecule 1 (ceacam1) (Figure 10). Most of these genes are involved with inflammation, complement activation, blood coagulation, and acute phage responses.

qPCR Verification Analysis
The gene expression relationship between the log2 of the RQ values from the RT-qPCR and the log2 of the transcript per million reads (TPM) values from the RNA-Seq was determined for 14 selected genes, including complement component c6 (c6), cxc chemokine receptor type 3 (cxcr3), galectin-3-binding protein a (igals3bp), glutathione s-transferase alpha tandem duplicate 1 (gsta4.1), hepcidine (hamp), interleukin 1 receptoe 2 (il1r2), interleukin 8 (cxcl8b/il8), bactericidal permeability-increasing protein (bpifcl), pentraxin-related protein ptx3 (ptx3a), ras-related protein orab-1 (orab1), amyloid protein a (app), suppressor of cytokine signaling 3a (socs3a), tumor necrosis factor receptor superfamily member 9 (tnfrsf9), and toll-like receptor 5 (tlr5a). As shown in Figure 11, there was a significant positive correlation between the RNA-Seq and RT-qPCR data. Correlation r 2 values ranged from 0.6 to 0.93 for all 14 genes. These results indicate that all examined DEGs were in agreement with the reference genome-guided RNA-Seq analysis. On the other hand, the qRT-PCR results for c6, bpifcl, igals3bp, and orab1 were not in good agreement with the de novo RNA-Seq analysis (S9 Figure). However, the qRT-PCR results of all other DEGs evaluated agreed with the de novo RNA-Seq analysis, with correlation R 2 values ranging from 0.6 to 0.95 for the other 10 genes ( Figure S9). Furthermore, we observed significant downregulation of genes encoding MHCII and IgM in all analyzed organs (Table 4). In addition, cd79a, cd79b, and cd209 were downregulated in the head-kidney, spleen, and liver (Table 4).

qPCR Verification Analysis
The gene expression relationship between the log 2 of the RQ values from the RT-qPCR and the log 2 of the transcript per million reads (TPM) values from the RNA-Seq was determined for 14 selected genes, including complement component c6 (c6), cxc chemokine receptor type 3 (cxcr3), galectin-3-binding protein a (igals3bp), glutathione s-transferase alpha tandem duplicate 1 (gsta4.1), hepcidine (hamp), interleukin 1 receptoe 2 (il1r2), interleukin 8 (cxcl8b/il8), bactericidal permeability-increasing protein (bpifcl), pentraxin-related protein ptx3 (ptx3a), rasrelated protein orab-1 (orab1), amyloid protein a (app), suppressor of cytokine signaling 3a (socs3a), tumor necrosis factor receptor superfamily member 9 (tnfrsf9), and toll-like receptor 5 (tlr5a). As shown in Figure 11, there was a significant positive correlation between the RNA-Seq and RT-qPCR data. Correlation r 2 values ranged from 0.6 to 0.93 for all 14 genes. These results indicate that all examined DEGs were in agreement with the reference genome-guided RNA-Seq analysis. On the other hand, the qRT-PCR results for c6, bpifcl, igals3bp, and orab1 were not in good agreement with the de novo RNA-Seq analysis ( Figure S9). However, the qRT-PCR results of all other DEGs evaluated agreed with the de novo RNA-Seq analysis, with correlation R 2 values ranging from 0.6 to 0.95 for the other 10 genes ( Figure S9).
Microorganisms 2022, 10, x FOR PEER REVIEW 21 of 33 Figure 11. Gene expression correlation between RT-qPCR and RNA-Seq data of 14 selected DEGs. RNA-Seq data are presented as log2TPM (X-axis), and RT-qPCR data are represented as log2RQ (Y-axis). Three different colours represent gene expression in the head-kidney (brown), spleen (red) and liver (green). The circle represents control samples; the square represents 3 dpi, and the triangle represents 10 dpi. Each symbol is an average of three fish at a particular time point in that tissue. RNA-Seq data are presented as log 2 TPM (X-axis), and RT-qPCR data are represented as log 2 RQ (Y-axis). Three different colours represent gene expression in the head-kidney (brown), spleen (red) and liver (green). The circle represents control samples; the square represents 3 dpi, and the triangle represents 10 dpi. Each symbol is an average of three fish at a particular time point in that tissue.

Discussion
Lumpfish is an emergent cleaner fish species in the North Atlantic region. However, diseases, including bacterial diseases, are affecting the performance of lumpfish and its extended utilization. A. salmonicida is a globally distributed pathogen that infects and kills lumpfish [3]. The infection kinetics of A. salmonicida in lumpfish and its response to early and lethal infection has not been described. In this study, we established a reproducible A. salmonicida systemic infection model in lumpfish. Additionally, we examined the transcriptome profile of internal organs, including the head-kidney, spleen, and liver of lumpfish injected with a lethal dose (10 4 CFU dose −1 ) of A. salmonicida, during early (3 dpi) and late infection stages (10 dpi). Head-kidney is known as a primary lymphoid organ as it is a hematopoietic tissue in the teleost, similar to the bone marrow of higher vertebrates [73]. B cell development, antigen-sampling and antigen retention have been described in teleost head-kidney [24,73,74]. The spleen is the primordial secondary lymphoid organ that contains macrophages, MHC class II+ cells, and T cells [24,73,75,76]. The liver is also an important organ that takes part in metabolism and defense [25,26,28,29], and it is also considered as a lymphoid organ, as non-parenchymal cells of the liver take part in antigen presentation and immunomodulatory functions. In addition, the liver encompasses large populations of natural killer cells and T cells [26,28]. This study analyzed the transcriptome response of the three main lymphoid tissues (head-kidney, spleen, and liver) of lumpfish during a lethal A. salmonicida infection.
The virulence of different A. salmonicida isolates varies in different fish hosts. For instance, A. salmonicida DH170821-10 showed relatively lower pathogenicity with an LD 50 of 6.4 × 10 3 CFU dose −1 in rainbow trout and coho salmon (Oncorhynchus kisutch) [77]. Another study described two highly pathogenic strains of A. salmonicida, MT1057, and MT423, with an LD 50 of 10 2 CFU dose −1 in Atlantic salmon but a lower virulence in halibut, with an LD 50 of 10 6 CFU dose −1 [78]. Our study showed that A. salmonicida J223 (Santander lab collection) is highly virulent for Newfoundland lumpfish. We determined that an ip infection of 10 2 bacterial cells per dose can kill 50 percent of the infected lumpfish population, which is similar to rainbow trout (Oncorhynchus mykiss), Chinese perch (Siniperca chuatsi), and Atlantic salmon [78,79]. The hyper-virulence of A. salmonicida J223 strain in lumpfish was further verified by another study conducted by our group, where a bath infection of lumpfish with 10 6 CFU mL −1 of A. salmonicida J223 caused 100% lethality within 14 dpi (unpublished data).
Subsequently, A. salmonicida infection kinetics in different organs was determined for different doses used to infect lumpfish. All lethal doses (10 3-10 5 CFU dose −1 ) showed the presence of A. salmonicida in the head-kidney at 3 dpi, suggesting that this organ is the primary A. salmonicida target, and from then it spreads to the spleen and liver, and finally, after 7 dpi, A. salmonicida infection in lumpfish becomes systemic (Figure 1). Similar to our findings, previous studies indicated that 3 to 4 days is a typical incubation period for A. salmonicida, where the bacterium rapidly disperses in the kidneys, followed by the spleen and liver [80,81]. Lumpfish infected with a low dose of A. salmonicida (10 1 CFU dose −1 ) established a persistent infection, as bacterial colonies were still detected after 30 dpi without causing mortalities. While A. salmonicida J223 strain lethal doses cause acute infections, in low doses it might cause chronic infections.
Similarly, Pseudomonas aeruginosa can cause both symptomatic acute and chronic infections. While acute infections often spread rapidly and can damage tissues as well as contribute to high mortality by sepsis, chronic infections can be carried on for years [82]. We did not explore the further mechanism of A. salmonicida mediated chronic infection here. Future studies to consider how A. salmonicida can utilize strategies to evade immune clearance to cause chronic infections would be helpful to explore the pathogenesis in marine teleosts.
To understand the transcriptome dynamics and their impact on gene expression levels, high-throughput RNA-Seq technology was used. RNA-Seq can effectively analyze transcript sequences and estimate gene expression levels that can be applied to the identification of DETs or DEGs between different experimental conditions [83,84]. RNA-Seq results in millions of short reads which need to be assembled into transcript sequences [85]. An RNA-Seq analysis allows for the distinguishing between individual transcripts (isoform) of a gene [85]. Analysis of DETs is essential in identifying differences between tissues [84]. However, the alignment of RNA-Seq reads to a certain gene allows researchers to study gene expression [86,87]. Gene expression estimation from the expression levels of transcripts provides more robust results [88]. Gene expression estimation allows researchers to determine DEGs under different conditions. Analyzing DEGs is more applicable for biological analysis, e.g., GO enrichment analysis [89]. Our study utilized two different approaches: de novo and reference-based, to assemble the transcriptome. With the availability of the reference genome, a reference-based assembly is more effective than a de novo assembly [90]; however, studies showed that the de novo assemblies were able to identify a complete gene content [55,60,[91][92][93][94][95]. We applied the de novo assembly approach at the isoform level that allowed us to determine DETs in 3 and 10 dpi of the head-kidney, spleen, and liver. However, the reference-based assembly approach allowed us to generate both DEGs and DETs using CLCGWB. Our results demonstrate that the total number of DETs identified by the de novo transcriptome assembly analysis was higher than the total number of DETs identified by the reference genome-guided transcriptome assembly analysis (5265 vs 4261, log 2 FC ≥ 1, FDR ≤ 0.05), which is similar to Kovi et al. [91] (Files S3 and S2). Intrinsic methodological issues of de novo analysis could generate misassembled transcripts [96]. The trinity de novo assembler might yield more transcripts because of lacking strand-specific information [97]. Subsequently, a BLAST search of all de novo DETs against the lumpfish genome identified that only 4839 (91.9%) de novo DETs are protein-coding transcripts. Hereafter, the corresponding genes of these de novo DETs were compared with the reference-based DEGs. We observed that only 25.9% of the genes were shared between the de novo and reference-based analysis, 17.3% of genes were unique in de novo analysis, and 56.8% genes were unique in reference-based analysis ( Figure 5).
In addition to this, our qPCR verification analysis demonstrates that the overall gene expression levels were underestimated by de novo analysis (Figures 11 and S9). Previous studies have found that the reference-based method surpasses the de novo method for characterizing transcriptome and gene expression [96,98,99]. Still, our study suggested that each method captured unique transcripts. Therefore, we adopted an integrative approach for GO enrichment analysis to bring more benefits for the better exploration of pathogenesis.
The number of DEGs and DETs was highest in the spleen, followed by the headkidney and liver at 3 dpi. Similarly, the number of DEG and DET were highest in the spleen, followed by liver and head-kidney at 10 dpi. However, in most cases, the number of DETs were lower than that of DEGs. This is because the gene-level expression is global. One gene can have several transcripts as a result of alternative splicing in eukaryotes and not all the DETs were significant (log 2 FC ≥ 1, FDR ≤ 0.05). Therefore, we cannot compare between the gene and transcript expression. Thus, moving forward, we used the gene-level analysis.
The head-kidney plays a key role in initiating the immune response in fish [73,100]. We observed that the initial inflammatory response was triggered in the head-kidney at 3 dpi (File S8). The histopathological analysis also detected the highest level of inflammation in the head-kidney, followed by the spleen and liver, respectively ( Figure 1F). Such responses correlate with the infection kinetics of A. salmonicida ( Figure 1B-E). Nevertheless, the spleen also was infected very fast and showed a tremendous amount of DEGs and enriched GO terms ( Figures 1C, 4D and 6C and Table 1). The spleen has a key role in promoting humoral immunity [101,102] and plays a key role in identifying cell damages [103]. This could be a reason for having the highest spleen response during A. salmonicida infection. The liver controls biochemical processes, including metabolism [104]. A metabolic arrest was suggested at 10 dpi in the liver (Files S6, S7 and S8). Interestingly, we observed fish lethargy (e.g., lack of appetite and swimming) starting at 7 dpi and continuing until death, which might relate to the metabolic arrest at the deadly point of the infection.
However, under some circumstances, these innate immune responses cause tissue damage and organ failure, eventually leading to death, which is a hallmark of sepsis [113]. During sepsis, the association of pattern recognition molecules with the pro-inflammatory mediators and activation of the NF-κB signaling cascade could cause the increased expression of proinflammatory cytokines [111]. Pro-inflammatory cytokines and complement components activate the coagulation cascade [114]. The coagulation system acts as a general host defense system to restrict the dissemination of pathogens by recruiting leukocytes, while fibrin promotes the adherence and migration of cells [115]. However, overactivation of the coagulation system during acute bacteremia causes disseminated intravascular coagulation (DIC), microvascular thrombosis-induced hypoxia, and a prolonged suppression of fibrinolysis, which contributes to multiorgan failure, abnormalities in host metabolism, immune suppression, septic like shock, and death [111,[115][116][117][118]. Interestingly, the downregulation of genes encoding hemoglobin subunits alpha and beta was identified in the spleen at 3 dpi, and the moribund lumpfish was visually noticed with the symptoms of hypoxia (Files S1, S7, and S8).
The unrestricted activation of inflammation, blood coagulation, and complement systems break the blood/tissue barrier and damage the host tissue and organs [119]. Interestingly, excessive hemorrhages in the lumpfish body, eyes, gills, or at the base of the fins, muscles, and organ tissues and astics were visually observed in moribund fish ( Figure S10). These observations suggest that detrimental and uncontrolled inflammation, overactivation of blood coagulation, and complement components might lead to a septiclike shock, which plays a significant role in the A. salmonicida mediated lethal infection of lumpfish. However, the septic response is an extremely complex reaction of inflammatory, anti-inflammatory, humoral, and cellular processes, and circulatory abnormalities, which are highly variable with the non-specific nature of the symptoms [120]. Therefore, further investigation is required to confirm sepsis in lumpfish.
Furthermore, il10 was upregulated in all three lumpfish organs at 10 dpi ( Figure 7). Previously, it was described that A. salmonicida elicits a significant increase in il10 expression in head-kidney leucocytes [121]. A similar effect was also described in Arctic char [81]. IL-10 can contribute to the immune suppression by inducing a Treg-mediated response. Deleting the T3SS genes of A. salmonicida decreases the host cytokine expression significantly [121]. A fully virulent A. salmonicida downregulated specific innate and adaptive immune gene expression and reduced the survival of the infected rainbow trout [122][123][124]. Consequently, we observed the downregulation of genes encoding MHCII and IgM in all analyzed organs (Table 4). Previous research on A. salmonicida infection in trout showed the downregulation of immunoglobulin light chains, constant and variable domains [108,109]. In addition, cd79a, cd79b, and cd209 were downregulated in the head-kidney, spleen, and liver (Table 4). CD79a and CD79b are B-cell antigen receptor complex-associated proteins α, and β chains play a crucial role in B cell development and antibody production [125]. CD209 is a C-type lectin, an essential PRR that participates in immune defense and microbial pathogenesis in mammals, and it is present on the surface of macrophages [126]. Coincident, previous studies on A. salmonicida infection in rainbow trout showed that cd209 was downregulated [108,109]. All of these observations suggest the A. salmonicida mediated immune suppression in lumpfish.
A. salmonicida virulence factor AopO is an ortholog of the Yersinia YopO/YpkA serine/theonine kinase. This serine/threonine kinase inhibits phagocytosis by blocking the Rac-dependent Fcγ receptor internalization pathway [122]. We observed the downregulation of the low-affinity immunoglobulin gamma Fc region receptor II-b-like (LOC117747925) in the spleen at 3 dpi (File S1), suggesting that A. salmonicida J223 might cause an antiphagocytic effect in lumpfish. However, this could be the result of an effect of undetermined signaling cascade, which needs further verification.
At least six A. salmonicida type-3 secretion system-related virulence factors, AexT, Ati2, AopH, AopO, AopN, and AopS could be responsible for disrupting the host cytoskeleton structure, which allows this pathogen to colonize and survive inside the host [30,122]. The GTPase activating domain and the ADP-ribosylating domain of AexT act on small monomeric GTPases of the Rho family (Rho, Rac, and Cdc42) and actin, respectively, and causes actin depolymerization and cell rounding [30,122]. Ati2 of Vibrio parahaemolyticus is responsible for the local detachment of the actin-binding proteins from the plasma membrane and induces membrane blebbing and cytolysis by hydrolyzing the host phosphatidylinositol 4,5-bisphosphate [30,122]. AopH, an ortholog of Yersinia YopH, is responsible for altering the actin cytoskeleton by dephosphorylating the tyrosine residue [30,122]. AopO, an ortholog of Yersinia YopO, prevents the actin distribution in the host cell [30,122]. AopN, an A. salmonicida effector, binds and sequesters αβ-tubulin and inhibits microtubule polymerization that induces mitotic arrest [30,122]. AopS, an ortholog of V. parahaemolyticus VopS, could inhibit the actin assembly by preventing the interaction of Rho GTPases with its downstream effectors [30,122]. In this study, lumpfish DEGs associated with cytoskeleton organization (e.g. actin-binding LIM protein 1-like, cdc42 effector protein 1b, rho GTPase-activating protein 4-like, rho guanine nucleotide exchange factor 10-like protein, rho guanine nucleotide exchange factor 18, rho-related GTP-binding protein RhoH, tubulin beta chain, tubulin monoglycylase TTLL3-like, tubulin polyglutamylase TTLL7, rasgrf2a, tppp3) were downregulated in the spleen at 10 dpi ( Figure 9 and Files S1, S7, and S8). In addition, downregulation of genes related to microtubule bundle formation was observed in the head-kidney at 10 dpi (e.g., genes encoding dynein assembly factors, dynein heavy chains, tppp3) Figure 9 and Files S1, S7, and S8). These findings indicate that the disruption of the lumpfish cytoskeleton might be possible by actin and microtubule depolymerization and mitotic arrest during A. salmonicida infection.
The A. salmonicida type-3 secretion system (T3SS) effector AopP induces apoptosis in affected cells by interfering with critical signal transduction pathways (i.e., NFκB signaling) that activate caspase 3 [30,122]. AopP hinders the NF-κB signaling pathway by restraining the transportation of the p50/p65 protein complex (NFKB1/RelA) into the target cell's nucleus [30,122], resulting in the septicemia and furuncles formation (subcutaneous wounds) in host tissue [127]. We observed the upregulation of several genes that positively regulate the apoptosis process in the head-kidney and spleen at 10 dpi (Files S1, S7, and S8). However, no caspases were differentially expressed in this study. We detected the upregulation of several regulators of the NFκB pathway, including κB (IκB) inhibitor, bcl-3, tnfrsf11b, aen, ddit4, nfkbia, and nuclear factor interleukin-3-regulated protein (nfil3) (Files S1 and S8). We did not observe the formation of furuncles in lumpfish skin that might be concurrent with no expression of caspases involved in apoptosis. A dual transcriptomic study and future in vitro experiments to identify the dysregulation of aopP of A. salmonicida in lumpfish lymphoid organs could be valuable for future research.
Certain bacterial pathogens could cause chronic inflammation and/or produce genotoxins that can damage proteins, lipids, metabolites, DNA, and RNA. For example, Helicobacter pylori infection downregulates DNA mismatch repair and base excision repair mechanisms [128]. The bacterial toxin can be a source of DNA double-strand breaks (DBSs), causing cell death [129]. DSBs induces DNA damage response (DDR), resulting in cell cycle arrest [130]. Our results indicate that A. salmonicida infection downregulates several genes involved in DNA damage/repair in the liver at 10 dpi (e.g., bard1, dna2, ercc1, ercc4, herc2, msh2, parp1, rad52, rad54l) ( Figure 10 and Files S1, S7, and S8). Therefore, several biological processes such as DNA replication, DNA and RNA metabolic processes, double-strand break repair, DNA repair, RNA metabolic process, gene expression, and cell cycle processes were enriched by the downregulated genes in liver 10 dpi (File S7). These findings suggests that A. salmonicida infection might provoke lumpfish DNA damage and cause cell cycle arrest in lumpfish liver.
Suitable biomarkers of sepsis and infection are necessary for monitoring fish disease conditions [131]. hamp, hp, app, ptx3, mmp13b, il1b, il8, il10, and il6 were significantly upregulated in the head-kidney, spleen, and liver of infected lumpfish, suggesting they could be used as biomarkers for the molecular diagnosis of A. salmonicida infection (Figure 7). Actually, most of these genes were suggested as biomarkers of sepsis in humans [132,133], suggesting a conserved response to septic shock in vertebrates. Genes encoding ras-related GTPase 1Ab, rho GTPase-related proteins, and microtubule-associated proteins might be proposed as biomarkers to identify A. salmonicida specific infection (Figure 9 and File S9). In addition, tlr5, c6, c7, fgb, fgg, f3a, socs3a, adora3, ceacam1, tppp3, tuba1a, ddit4, dna2, and msh2 can also potentially be proposed as a biomarker to detect A. salmonicida lethal infection in lumpfish. Multiplex qPCR assays for these genes could then be developed to detect early A. salmonicida infection in lumpfish. These high-throughput technologies could accelerate the identification of potential biomarkers for various diagnostic and therapeutic developments in the future lumpfish aquaculture industry and explore the repose to septic shock in marine teleost.

Conclusions
A. salmonicida has evolved a myriad of mechanisms to counteract and modulate the host responses. Only 10 2 cells of A. salmonicida can kill 50% of the lumpfish population. Overall, our study characterizes A. salmonicida infection kinetics in lumpfish head-kidney, spleen, and liver ( Figure 1) and proposes an infection model for lumpfish molecular responses at the early and lethal point of infection (Figure 12). The model suggests that A. salmonicida might induce lethal infection in lumpfish by uncontrolled and detrimental blood coagulation, complement activation, and inflammation. Such responses could lead to hypoxia, internal organ hemorrhages, suppression of the adaptive immune system, and impairment of the DNA repair system, which results in cell cycle arrest, and, ultimately, death ( Figure 12). Also, A. salmonicida might destabilize the cytoskeleton structure by depolymerizing actin and microtubules to colonize and survive inside the lumpfish ( Figure 12). In addition, A. salmonicida may be capable of inhibiting the NF-kB signaling pathway and the induction of the apoptotic process ( Figure 12). These findings could help a global understanding of the molecular network of A. salmonicida-lumpfish host interactions, which is essential for developing effective treatments. Furthermore, our analysis provides a guideline for future experimental designs to study A. salmonicida pathogenesis in lumpfish.
Microorganisms 2022, 10, x FOR PEER REVIEW 27 of 33 depolymerizing actin and microtubules to colonize and survive inside the lumpfish (Figure 12). In addition, A. salmonicida may be capable of inhibiting the NF-kB signaling pathway and the induction of the apoptotic process ( Figure 12). These findings could help a global understanding of the molecular network of A. salmonicida-lumpfish host interactions, which is essential for developing effective treatments. Furthermore, our analysis provides a guideline for future experimental designs to study A. salmonicida pathogenesis in lumpfish. Figure 12. Aeromonas salmonicida infection model in lumpfish lymphoid organs. A. salmonicida early (3 dpi) and late (10 dpi) infection model suggests that A. salmonicida induces lethal infection in lumpfish by uncontrolled and detrimental blood coagulation, complement activation, and inflammation. Such responses lead to hypoxia, internal organ hemorrhages, suppression of the adaptive immune system, and impairment of the DNA repair system, which results in cell cycle arrest, and, ultimately, death. Furthermore, A. salmonicida could destabilizes the cytoskeleton structure by depolymerizing actin and microtubule.