Forty-eight samples of head kidney from three time points (0 dpi, 517 ddpv; 8 dpi, 605 ddpv; 29 dpi, 836 ddpv) consisting of three experimental conditions (Sham control, FM-vaccinates, FM+R-vaccinates; n = 5–7 individuals per condition) were selected for library construction, RNA-seq analysis, and RT-qPCR analysis.

RNA was extracted from frozen head kidney samples using Tri-Reagent as previously described (26), and genomic DNA contamination was eliminated after DNase treatment (Ambion® TURBO DNA-free™). Quality of resulting purified RNA was determined using Experion™ Automated Electrophoresis Station (BioRad) and Nanodrop 3,000 (Thermo Fisher) was used to test both purity and quantity. RNA was stored at −80°C until subsequent use.

DNA was isolated from 81 individual head kidney tissues (3 fish per group pre-challenge, and 4 fish per tank [triplicate tanks × 3 treatment groups × 2 post-challenge samplings]) using a Qiagen® DNeasy Blood and Tissue Kit following manufacturer's instructions with the following changes: samples (~25 mg) were homogenized for 20 min at 50 Hz (TissueLyser, Qiagen) and incubated at 56°C for 4 h. Resulting purified DNA was quantified spectrophotometrically (Nanodrop 3000) and diluted with nuclease-free water to 50 ng/ml. Integrity of gDNA was assessed on a 1% agarose gel before downstream enzymatic assays.

Posterior kidneys were cultured from fish collected on designated sampling days and mortalities occurring throughout the trial by aseptically plating culture swabs on BA media at 13°C. The plates were checked at least once weekly for bacterial growth. The presence of bacterial growth, and the characteristic brown colonies on tyrosine-rich media was considered a positive criterion for typical Asal. Plates showing uncharacteristic growth (possible contamination) or no growth were considered negative for the presence of the bacteria.

Serial sections of paraffin-embedded formalin-fixed tissues were sectioned at 7 μm and either stained with haematoxylin and eosin (H&E) or Wright's Giemsa following standard histological techniques or probed with anti-trout IgM monoclonal antibodies using the Expose Mouse-Rabbit Specific HRP/DAB kit (Abcam) following manufacturer's instructions. Briefly, after deparaffinization and rehydration, sections were heated to 100°C in antigen retrieval buffer (pH 9) cooled to room temperature for 10 min in phosphate-buffered saline (PBS), and then washed twice in Tris-buffered saline plus 0.2% Tween 20 (TBS-T; pH 8.0) for 5 min with gentle agitation. Sections were blocked in protein blocker for 10 min before gentle rinsing with TBS-T. The sections were incubated with α-trout mIgM (1:50) in TBS-T and 1% bovine serum albumin (BSA; Sigma) overnight at 4°C in a humid chamber. After incubation, the sections were washed in TBS-T (two times for 5 min each) and incubated in a mouse-specifying reagent (Expose mouse/rabbit specific horseradish peroxidase [HRP]/3,3 = -diaminobenzidine [DAB] kit; Abcam) for 10 min, followed by a 10-min incubation in hydrogen peroxide blocker. Labeled cells were detected after a 15-min incubation in a goat anti-rabbit HRP conjugate followed by 10 min with DAB in PBS with 0.015% H₂O₂. All sections were counterstained in 1% Alcian blue (3 min) and Mayer's hematoxylin (diluted 1/20, 30 s), dehydrated in graded ethanol, cleared in xylene, and mounted (Permount). Sections treated with irrelevant antibodies served as negative controls, while sections known to contain IgM-labeled cells served as positive controls. After staining or immunolabeling, sections were visualized with a ZEISS AxioCam IC.

Total RNA was submitted for library construction and sequencing services by the McGill University and Génome Québec Innovation Center, Montréal, Canada.

RNA quality and purity were assessed using an Experion Bioanalyzer (BioRad), and only samples with a minimum RIN of 6 proceeded to library construction. Forty-eight TRUseq 100 bp pair-end stranded mRNA libraries were generated from the experimental conditions listed in Supplemental Table 1 and were sequenced in 4 lanes (12 samples per lane) on an Illumina HiSeq 2000 platform. De novo transcriptome assembly was performed on resulting reads following the pipeline described by Haas et al. (27) based on the Trinity assembly software suite v 2.1 (28). Briefly, reads were trimmed with Trimmomatic software (29) with a minimal Phred score of 33 and a minimal length of 32 bp. A normalized metric of reads was generated using Trinity normalization utility and surviving paired reads were assembled using the Trinity assembler (27). The final assembly quality was checked using R-correct and contigs with poor read support were removed. Trinotate v2.0.2 was used to identify putative coding transcripts and all putative transcripts were aligned against the UniProt protein database (downloaded March 2016) using the blastx program from the NCBI BLAST family as implemented in the blast+ package. Annotation was assigned to each longest putative coding transcript derived from each de Bruijn graph component (unigene) based on highest blastx score with an Expected value (E) cut-off <1e-10.

Read alignment to the reference transcriptome was performed using RNA-Seq by Expectation Maximization Method (RSEM) described by Li et al. (30) using the de novo assembled transcriptome as a reference. This method was selected as it allows for the inclusion of non-uniquely mapping reads in abundance estimates through the use of a likelihood method that accounts for read mapping uncertainty.

RSEM quantities were used as the input for the programs edgeR (31) and DESeq2 (32) using scripts provided with the Trinity pipeline and a list of DE transcripts representing the mean expression values from samples for every treatment and time combination. Lists of DE transcripts generated using DESeq2 were generally larger and were almost always contained within the edgeR DE lists, and for this reason we chose to use the more conservative edgeR-generated DE lists going forward with analysis. Differential expression was tested at a significance level of α = 0.05. Transcripts were considered differentially expressed at a p < 0.05, and log₂fold-change > 2, following a false discovery rate (FDR) adjustment of 5% (0.05) as implemented by the BH algorithm in the p.adjust function of the R stats package.

Functional gene ontology was conducted using GOrilla (33) by comparing groups of differentially expressed genes to the de novo reference assembly as the background with an FDR-corrected p < 0.001. Using this open-access software, enriched GO terms are identified from ranked lists of DEGs as well as from target lists of DEGs compared to a background. Resulting enriched lists were summarized and visualized using REViGO (34) and Cytoscape.

To validate sequencing data, gene-specific primers were designed using Primer3 (35) and the contig sequences for nine significantly expressed transcripts, and analyzed through qPCR (Supplemental Table 2). Total RNA was extracted from samples of head kidney as described above. For every sample, 750 ng of DNase-treated total RNA was reverse transcribed using a High Capacity cDNA Reverse Transcription Kit (Thermo Fisher) following manufacturer's instructions in a reaction volume of 20 μl. qPCR was performed in triplicate reactions containing 2X Universal SYBR Green Master Mix (BioRad), 100 nm of each forward and reverse primer, and 2 μl template (diluted 1/20) in a total reaction volume of 11 μl. Cycling was performed on a CFX96 thermal cycler (BioRad) with the following thermal profile: 95°C for 30 s (1 cycle), 95°C for 15 s then 60°C for 30 s (40 cycles), followed by melt curve analyses from 65 to 95°C with fluorescence being read every 0.5 s with a ramp rate of 0.5°C to ensure amplification specificity. An equal volume of every cDNA sample was pooled and diluted 10-fold in five-steps to determine amplification efficiency and linearity for every primer pair. To confirm an absence of gDNA contamination, noRTs (no reverse transcriptase control) were performed for each gene.

Target gene expression was normalized to the two most stable of nine reference genes: 60S ribosomal protein L7 (60S), beta-2-microglobulin (β2m), elongation factor 1-alpha (efIα), glyceraldehyde-3-phosphate dehydrogenase (gapdh), hypoxanthine-guanine phosphoribosyltransferase (hprt), ribosomal protein S20 (rps20), beta-actin (ß-actin), eukaryotic initiation factor 5A(eif5), and tubulin alpha chain (tubulin-α). Elongation factor 1-alpha and rps20 showed the highest stability, with geNorm M value and coefficient of variation of 0.993 and 0.37, respectively, using the geNorm algorithm found in qBASE⁺ software (36). Calibrated normalized relative quantities (CNRQ) were calculated with sample-specific normalization factors and gene-specific amplification efficiencies, and internal positive controls were included to calibrate run-to-run variation among plates.

Amplicons for each gene were gel purified using QIAquick Gel Extraction Kit (Qiagen) as per the manufacturer's instructions. Purified products were sent to The Center for Applied Genomics (TCAG) at the Hospital for Sick Children (Toronto, Canada) for primer verification (Supplemental Table 3).

Specific primers for the Asal A449 aopO gene located in the low-copy-number pAsa5 plasmid (GenBank accession no. DQ386862.1) have been published previously (37). Extracted DNA from posterior kidneys of moribund and sampled fish was used to confirm Asal infection by conventional and real-time PCR. PCR products were synthesized using 12.5 μl GoTaq (ProMega), 2 μl of forward aopO primer (5′-AGCTCATCCAATGTTCGGTATT-3′), 2 μl of reverse aopO primer (5′-AAGTTCATCGTGCTGTTCCA-3′), 5 μl of DNA template (50 ng/μl) and 3.5 μl of nuclease-free water. The samples were run on an Eppendorf MasterCycler thermal cycler with the following thermal profile: initial heating at 95°C for 2 min, denaturation at 95°C for 1 min, annealing at 60°C for 30 s, extension at 72°C for 30 s, repeat for 34 cycles, final extension at 72°C for 5 min, and hold at 4°C. DNA from pure cultures of Asal was included as a positive control, whereas reactions lacking DNA template were used as negative controls. To confirm if samples were positive or negative for the aopO gene, PCR products were separated on a 2% agarose gel, and visualized using an ultraviolet light transilluminator (UV Transilluminator 2000, Bio Rad), with the presence of a 119-bp band considered as a positive result.

Real-time PCR was performed on a CFX96 thermal cycler (BioRad) and the CFX Manager™ software detection system (BioRad). The reaction mix contained 12.5 μl of 2X Universal SYBR Green Master Mix, 400 nM of each primer, 2 μl of template, and nuclease-free water to reach a final volume of 25 μl. The thermal cycling profile was comprised of an initial incubation of 95°C for 30 s, 35 cycles of 95°C for 15 s then 62°C for 30 s. Melt curve analysis was performed to ensure amplification specificity using a temperature gradient from 65 to 95°C with fluorescence being read every 0.5 s with a ramp rate of 0.5°C/s. Confirmed positive template were used as a positive control, while DNA extracted from fish prior to Asal-challenge were used as negative controls. All real-time PCR assays were performed in triplicate, and an NTC was included on every plate.

The Kaplan-Meier (KM) survival curves were calculated in R using the survival package (R Core Developmental Team, V3.4.2).

The relative proportion survival (RPS) was calculated for the average of the replicates for each combination of treatment factors following the formula (38):

where a = cumulative mortalities in vaccinates

b = cumulative mortalities in non-vaccinates

ELISA data were analyzed using a two-way ANOVA followed by Tukey HSD. Correlations between mortality and antibody titer were detected using the Pearson method and performed in R (V3.4.2).

An average of 15.7 million reads per library was generated, with 73% of all reads mapping to the reference transcriptome (Supplemental Table 5). All contigs from the de novo assembly were analyzed using Transdecoder to identify the best putative ORF. Resulting protein sequences were annotated using the Trinotate v2.0.2 using the uniport database retrieved on March 2016.

RNA-seq libraries were generated from 48 samples of Arctic charr arbitrarily chosen from 9 discrete treatment conditions (Supplemental Table 1). Raw sequence data from all libraries has been deposited in the NCBI Sequence Read Archive as accession PRJNA507334.

Transcriptome de novo assembly yielded a total of 664,663 contigs with an average transcript length of 786 bp. These sequences represent 423,594 trinity genes, of which, ~181,321 (42.8%) were suggested as being protein-coding following Trinotate and Uniprot blastx annotation analysis (Supplemental Data File 1). Of these, 69,019 were annotated by UniProt. Abundance estimates (RSEM expected counts) for all transcripts identified in each RNA-seq library is provided in Supplemental Data File 2. To provide concordance in differential expression we first used two different pipelines (edgeR, DESeq2). Lists of DE transcripts generated from DESeq2 were generally larger and fully encompassed edgeR DE transcript lists (Supplemental Figure 2). We decided to consider data obtained by the more conservative edgeR pipeline for subsequent analysis. Lists of differentially expressed transcripts identified from abundance estimates using edgeR are provided in Supplemental Data File 3.

The transcriptomic profile of both groups of vaccinated Arctic charr followed a similar pattern to that of Sham controls, characterized by low numbers of differentially regulated transcripts prior to bacterial challenge, followed by an increase in over-expressed transcripts during the bacterial infection, ending with a decrease in over-expressed transcripts in post-challenge, surviving fish. However, in contrast to Sham controls, the degree and magnitude of the response in vaccinated fish was significantly reduced. Furthermore, the transcriptomic response was observed to differ depending on the vaccine administered with a much stronger transcriptomic response in FM-vaccinates compared to FM+R-vaccinates (Figure 6; Supplemental Tables 6, 7).

To evaluate the effect of vaccination on the transcriptome, we first compared the transcriptomic profiles of vaccinated (FM or FM+R) to unvaccinated (Sham) fish prior to infection with Asal (517 ddpv). There were 145 and 160 differential expressed genes (DEGs) in fish vaccinated with FM+R and FM, respectively, with 19 transcripts concordantly overexpressed by both vaccinate groups (Table 2). Three of these transcripts (40S ribosomal protein S12, ubiquitin carboxyl-terminal hydrolase 3, and 7SK snRNA methylphosphate capping enzyme) were expressed in opposing directions (i.e., up-regulated in one group while down-regulated in the other).

Compared to Sham controls, the transcriptomic response of FM-vaccinates at 517 ddpv was characterized mainly by over-expression of genes involved in the acute phase response (e.g., alpha-2-HS-glycoprotein, serum albumin 1, ladderlectin), complement and coagulation cascade (e.g., complement factor H-related protein 1, complement C1r subcomponent, complement C3, alpha-2-macroglobulin; Figure 7), and iron homeostasis (e.g., hemopexin, serotransferrin-2, ceruloplasmin). In contrast, FM+R-vaccinates were characterized by up-regulation of genes involved in protein synthesis (e.g., 60S ribosomal protein L13a, 40S ribosomal protein S12), cellular transport (e.g., choline transporter-like protein 4, phospholipid-transporting ATPase ID) and transcription (e.g., eukaryotic translation initiation factor 3 subunit L, nuclear receptor coactivator 4, CREB-regulated transcription coactivator 2).

There was a high number of over-expressed transcripts associated with complement and coagulation in the FM-vaccinated charr (n = 13), compared to the group receiving the combined vaccination of FM+R (n = 3; Supplemental Table 8). Additionally, we only detected significant up-regulation of transcripts involved in the acute-phase response in FM-vaccinates prior to infection with Asal, such as alpha-2-HS-glycoprotein (Fetuin-A, FC = 17.3), serum albumin 1 (FC = 13.6), ladderlectin (FC = 13.1), apolipoprotein B-100 (FC = 12.5), and apolipoprotein A-I-1 (FC = 12.4).

Compared to PBS-injected controls, at 8 dpi/605 ddpv, the number of differentially expressed transcripts in FM-vaccinated fish was greater than the response in FM+R-vaccinated fish (1,163 and 145, respectively; Figure 6C, Supplemental Table 7). Of these, only 52 were concordantly over-expressed in both groups, including tissue remodeling enzymes (matrix metalloproteinase-9), hemostasis-associated genes (heme-binding protein 2, hemoglobin subunit alpha), mitochondrial enzymes ([3-methyl-2-oxobutanoate dehydrogenase [lipoamide]] kinase), cytoskeleton-associated transcripts (keratin type I cytoskeletal 18, WAS/WASL-interacting protein family member 1), and actin-binding enzymes (gelsolin, dematin) (Supplemental Table 9).

Transcripts specifically upregulated in FM-vaccinated fish were characterized by enrichment in “RNA metabolic process” (GO:0007275), “regulation of gene expression” (GO:0010468), and “multicellular organism development” (GO:0035987). In contrast, in FM+R-vaccinated fish there was enrichment of “granulocyte chemotaxis” (GO:0071621), and “dynamin family protein polymerization involved in mitochondrial fission” (GO:0003374) (Figure 8; Supplemental Data File 4).

There was no enrichment in down-regulated transcripts in FM-vaccinated Arctic charr, however, in FM+R-vaccinates there was enrichment in “respiratory electron transport chain” (GO:0022904), and “regulation of transcription of nucleolar large rRNA by RNA polymerase I” (GO:1901836).

We observed over-expression of several adaptive immunity-associated transcripts in immunized fish at 8 dpi, with a higher number of these transcripts over-expressed in FM-vaccinates. These included BCL-6 corepressor-like protein 1 (FC = 8.54), Ig heavy chain V region (FC = 5.44), Ig heavy chain V region BCL1 (FC = 8.32), and Ig kappa chain V-I region DEE (FC = 7.44). One transcript associated with adaptive immunity (Ig lambda-3 chain C region) was concordantly over-expressed as multiple isoforms in both FM- and FM+R-vaccinated fish (FC = 6.52–8.36 and FC = 3.6–8.95, respectively).

Interestingly, compared to unvaccinated controls, there was significant upregulation of a hemolytic protein, neoverrucotoxin subunit alpha, only in the FM-vaccinated group (FC = 8.9–10.3).

At 29 dpi/836 ddpv, the number of over-expressed transcripts in the FM-vaccinates decreased to a comparable level to that of FM+R-vaccinates (129 and 90, respectively; Figure 6C, Supplemental Table 7). Fifteen transcripts were concordantly over-expressed between vaccinated groups, with functions in RNA-binding (RNA-binding protein 10), DNA-binding (coiled-coil and C2 domain-containing protein 1A), DNA-repair enzymes (tyrosyl-DNA phosphodiesterase 2), transferases (probable E3 ubiquitin-protein ligase HERC1), helicases (fanconi anemia group J protein homolog), nuclear protein transport (importin-7), microtubule organization (centrosomal protein of 170 kDa), signal transduction (leucine-rich repeat-containing protein 28), and cytoskeletal organization (bromodomain and WD repeat-containing protein 3) (Table 3). Within these genes, three were differentially expressed between vaccine groups (transposon TX1 uncharacterized 149 kDa protein, synergin gamma, and membralin).

Transcripts specifically up-regulated in FM-vaccinated fish included those involved in DNA repair (e.g., DNA cross-link repair 1A protein, fanconi anemia group J protein homolog, tyrosyl-DNA phosphodiesterase 2), protein phosphatase/kinase activity (e.g., protein phosphatase 1H, calmodulin), transcriptional regulation (e.g., coiled-coil and C2 domain-containing protein 1A, breast cancer metastasis-suppressor 1-like protein-A, polyhomeotic-like protein 1), cell cycle regulation (e.g., aurora kinase), and cellular trafficking (e.g., probable E3 ubiquitin-protein ligase HERC1, importin-7). GO analysis detected enrichment in processes such as “negative regulation of phospholipase activity” (GO:0010519), and “epithelial fluid transport” (GO:0042045) (Figure 8; Supplemental Data File 4).

The transcriptional profile of FM+R-vaccinates also included many genes involved in the same pathways; however, in addition, we observed significant up-regulation in cellular respiration (e.g., NADH-ubiquinone oxidoreductase 75 kDa subunit, mitochondrial, adenylosuccinate lyase, glucokinase, COX assembly mitochondrial protein 2 homolog), adaptive immune regulation (e.g., butyrophilin subfamily 1 member A1, nuclear factor of activated T-cells 5), and proinflammatory regulation (e.g., leukotriene A-4 hydrolase).

There was significant downregulation at 29 dpi in FM-vaccinates, with enrichment in processes such as “protein activation cascade” (GO:0072376), “acute-phase response” (GO:0006953), “platelet degranulation” (GO:0002576), and “complement activation” (GO:0006956) (Supplemental Data File 4). In contrast, down-regulated genes in FM+R-vaccinates were only enriched for one biological process, “cellular senescence” (GO:0090398).

There were 1,334 transcripts up- and 511 transcripts down-regulated in Asal-infected Sham controls at 8 dpi compared to uninfected Sham controls. Functional Gene Ontology (GO) of upregulated transcripts revealed enrichment of biological processes involved in a number of physiological functions including “cytokine-mediated signaling pathway” (GO:0019221), “regulation of cell proliferation” (GO:0042127), and “response to lipopolysaccharide” (GO:0071222) (Figure 8; Supplemental Data File 3). In contrast, enrichment in downregulated transcripts included processes such as “secretion” (GO:0046903), “exocytosis” (GO:0006887), and “regulation of reactive oxygen species metabolic process” (GO:20000377).

As the complement and coagulation cascade was highly upregulated in vaccinated fish prior to Asal-challenge, we were interested in comparing the expression profiles of these transcripts in unvaccinated fish over the experimental time period. Interestingly, there was high concordance in expression of many transcripts between non-vaccinated surviving charr and infected vaccinated charr (Table 4).

Compared to uninfected controls, the transcriptomic response decreased in surviving unvaccinated charr at 29 dpi, with 123 up-regulated transcripts and 130 down-regulated transcripts. Gene ontology analysis of up-regulated transcripts indicated enrichment in biological processes involved in cell-cycle process and proliferation, including “cell cycle phase transition” (GO:0044770), “mitotic cell cycle process” (GO:1903047), and “meiotic spindle organization” (GO:0000212) (Supplemental Data File 4). Enrichment of down-regulated transcripts resulted in one biological process of “negative regulation of erythrocyte differentiation” (GO:0045647) (Supplemental Data File 4).

We then compared the transcriptomic response between the two vaccine-groups in attempts to identify differences in the host response that might account for the increased protection in the FM+R-vaccinates.

Compared to Sham controls prior to bacterial challenge (517 ddpv), there were 125 over-expressed transcripts, with 47 up-regulated in FM-vaccinated charr compared to 78 up-regulated in the FM+R-vaccinated group. Gene Ontology analysis of up-regulated transcripts in FM-vaccinates revealed enrichment in biological processes such as “acute phase response” (GO:0006953), “regulation of proteolysis” (GO:0030162), and “platelet degranulation” (GO:0002576). In contrast, FM+R-vaccinates showed no enrichment of GO terms. Enrichment of downregulated transcripts at 517 ddpv in FM-vaccinates revealed enrichment in processes including “positive regulation of lamellipodium organization” (GO:1902743), and “cellular response to biotic stimulus” (GO:0071216). In FM+R-vaccinates there was one category enriched in downregulated transcripts, “negative regulation of chromatin binding” (GO:0035562).

During the proliferative growth stage of Asal infection (8 dpi, 605 ddpv), 108 transcripts were over-expressed, with 62 up-regulated in FM-vaccinated fish, compared to 46 up-regulated transcripts in the FM+R-vaccinated fish. Functional GO of upregulated transcripts revealed FM-vaccinates had enrichment in “defense response” (GO:0006952), “positive regulation of establishment of T cell polarity” (GO:1903905), and “inflammatory response” (GO:0006954). In contrast, FM+R-vaccinates had enrichment of “androgen receptor signaling pathway” (GO:0030521), “positive regulation of actin nucleation” (GO0051127) (Supplemental Data File 4). In contrast, GO enrichment of downregulated transcripts of FM-vaccinated fish included processes of “chemotaxis” (GO:0006935), “regulation of signaling” (GO:0023051), and “leukocyte migration” (GO:0050900). In FM+R-vaccinates, there was much fewer enriched categories in downregulated transcripts, including “positive regulation of phagocytosis” (GO:0050766), “positive regulation of endocytosis” (GO:0045807), and “PERK-mediated unfolded protein response” (GO:0036499).

The response in surviving vaccinated fish (29 dpi, 836 ddpv) reduced to 85 over-expressed transcripts comprising 39 up-regulated in FM-vaccinates, and 46 up-regulated transcripts in FM+R-vaccinates. In FM-vaccinated fish, there was enrichment in upregulated transcripts for processes like “histone H3-K9 modification” (GO:0006413:000061647), “regulation of catabolic process” (GO:0009894), and “peptidyl-lysine monomethylation” (GO:0018026). There was no enrichment of GO terms FM+R-vaccinates in upregulated transcripts, and there was no enrichment in downregulated transcripts for either FM- or FM+R-vaccine group.

Comparison between overexpressed transcripts during the experimental time course showed that there was a similar profile between 517 ddpv and 8 dpi (~86% similarity). There were 16 transcripts specific to the pre-challenge FM-vaccinated group with significant activation in several humoral immune effectors (Table 5). In contrast, the profile of overexpressed transcripts in surviving vaccinates was discrete.

Nine genes of interest (B-cell linker protein, mx2, il6rb, c7, ladd, fibronectin, fibrinogen gamma chain, hsp90b) were chosen for validation with RT-qPCR on the RNA-seq samples, as well as on an expanded number of samples to increase the power of testing. Some of these genes were chosen to cross validate the DE results from RNA-sequencing, and others were of immunological importance. Fold-change comparisons between both methods revealed a significant positive correlation (Pearson's r = 0.60 p = 0.024) (Supplemental Data File 5).

We then correlated the expression of these immune associated genes with the expression of apoO from Asal. There was significant negative correlation in both blnk and finc, indicating that presence of Asal negatively affected expression of these genes (Supplemental Table 10).

Acute phase response genes ladderlectin, complement component 3, and complement component 7 were expressed higher in FM-vaccinates than in PBS-injected controls as identified by RNA-seq (Figure 9A). Genes involved in immunological processes such as B-cell differentiation and cytokine signaling were also confirmed to be differentially regulated through vaccination. ForteMicro®-immunized fish had a significantly higher expression of finc and blnk compared to Sham controls after infection with Asal (8 dpi, 605 ddpv; Figure 9B), with expression increasing in the Sham and FM+R group by the final sampling time. Other genes examined using RT-qPCR, as mentioned above, had concordant expression with the RNA-seq data, but were not all significant over time and treatment comparisons (Supplemental Figure 3).