Comparative Transcriptome Analysis of Spleen Reveals Potential Regulation of Genes and Immune Pathways Following Administration of Aeromonas salmonicida subsp. masoucida Vaccine in Atlantic Salmon (Salmo salar).

Aeromonas salmonicida is a global fish pathogen. Aeromonas salmonicida subsp. masoucida (ASM) is classified as atypical A. salmonicida and caused huge losses to salmonid industry in China. Hence, it is of great significance to develop ASM vaccine and explore its protection mechanism in salmonids. In this regard, we conducted RNA-seq analysis with spleen tissue of Atlantic salmon after ASM vaccination to reveal genes, their expression patterns, and pathways involved in immune protections. In our results, a total of 441.63 million clean reads were obtained, and 389.37 million reads were mapped onto the Atlantic salmon reference genome. In addition, 1125, 2126, 1098, 820, and 1351 genes were significantly up-regulated, and 747, 2626, 818, 254, and 908 genes were significantly down-regulated post-ASM vaccination at 12 h, 24 h, 1 month, 2 months, and 3 months, respectively. Subsequent pathway analysis revealed that many differentially expressed genes (DEGs) following ASM vaccination were involved in cytokine-cytokine receptor interaction (TNFRSF11b, IL-17RA, CCR9, and CXCL11), HTLV-I infection (MR1 and HTLV-1), MAPK signaling pathway (MAPK, IL8, and TNF-α-1), PI3K-Akt signaling pathway (PIK3R3, THBS4, and COL2A1), and TNF signaling pathway (PTGS2, TNFRSF21-l, and CXCL10). Finally, the results of qRT-PCR showed a significant correlation with RNA-seq results, suggesting the reliability of RNA-seq for gene expression analysis. This study provided insights into regulation of gene expression and their involved pathways in Atlantic salmon spleen in responses to vaccine, and set the foundation for further study on the vaccine protective mechanism in Atlantic salmon as well as other teleost species.

Introduction

Atlantic salmon (Salmo salar), one of the most economically important fish utilized for aquaculture, is a key species in marine and freshwater ecosystems and has a high economic value to ecotourism and fisheries (Andrew et al. 2021; Wang et al. 2019). The economic importance of this species as well as other factors have led to it becoming a model species for biological research (Andrew et al. 2021). Atlantic salmon is by quantity the largest species of salmonids produced worldwide (FAO, 2018), and they are mainly farmed in cold waters in Northern Europe, North America, Chile, and Tasmania (Gjessing et al. 2020), and were introduced into China using recirculating aquaculture system (RAS) in 2010 due to their high nutritional and economic value (Wang et al. 2019). The anadromous salmon life cycle is mimicked in aquaculture, where juvenile fish are farmed in indoor fresh water facilities, artificially smoltified, and then transferred to sea water for the grow-out phase (Gjessing et al. 2020).

Aeromonas salmonicida is an important bacterial fish pathogen causing furunculosis in many different freshwater and marine fish species worldwide (Dallaire-Dufresne et al. 2014) and causes significant economic losses to global Atlantic salmon industry (Janda and Abbott 2010). A. salmonicida is a non-motile psychrophilic species made up of five subspecies (salmonicida, achromogenes, masoucida, pectinolytica, and smithia) (Garrity 2001). A. salmonicida subsp. salmonicida (ASS) was classified as “typical” and any isolate deviating phenotypically was regarded as “atypical” (Reith et al. 2008). In a recent study, A. salmonicida subspecies were divided into either of 14 distinct clusters or A-layer types based on partial nucleotide sequences of the virulence array protein (vapA) gene (Gulla et al. 2016). A. salmonicida subsp. masoucida (ASM) strains, isolated from China, Japan, and South Korea, were clustered in cluster VII among these 14 clusters. Several ASM strains were recently isolated from aquaculture farm of salmon in the city of Yantai, China, and were reported to bring huge economic damages to the local salmonid aquaculture (Du et al. 2015; Yan et al. 2021). Although the host response to the bacterium and vaccination has been well documented for A. salmonicida (Gudmundsdottir et al. 2010; Midtlyng et al. 1996a; Villumsen et al. 2012; Du et al. 2015), the effective ASM vaccine is still lacking, and little is known on the molecular mechanisms of host immune response after ASM vaccination.

Fish maintains a healthy state by defending itself against environmental pathogens by a complex system of innate defense mechanisms (Magnadóttir 2006). The spleen is regarded as the primordial secondary lymphoid organ, and almost all gnathostomes possess this organ in which adaptive immune responses are generated, concentrating the antigens for interactions between antigen-specific T cells, B cells and antigen-presenting cells (APCs) (Flajnik 2018; Bjørgen and Koppang 2021; Qi et al. 2016). Teleost spleen is the major source of immunoglobulin and play vital roles in preventing further lesions during bacterial infections (Li et al. 2015). In addition, the teleost spleen was reported to be capable of initiating the immune responses associated with calcium (Hwang et al. 2015), and producing immune globulin (Grove et al. 2006) and various lytic enzymes (Milla et al. 2010) to protect fish against pathogens.

Transcriptome, as an effective analysis tool, enables an in-depth study of complicated physiological pathways, including development, immune responses, physiology, cellular fate, and disease progression (Maekawa et al. 2019; Wang et al. 2009). Recently, transcriptomics technologies have been fully utilized to explore the participations of teleost spleen in response to different pathogens (Ali et al. 2014; Papetti et al. 2015; Kim et al. 2020). However, the information on the roles of Atlantic salmon spleen during immune activities at the transcriptional level is quite limited, even if the transcriptional responses were explored in Atlantic salmon in several tissues including blood, liver, head kidney, and gill (Andrew et al. 2020; Martin et al. 2006; Botwright et al. 2021). To expand the knowledge at the molecular level, we conducted RNA-seq analysis to identify the differentially expressed genes (DEGs) and the immune signaling pathway they involved in Atlantic salmon spleen after intraperitoneal injection of ASM vaccine. These findings will contribute to fish immunotherapy for the prevention and treatment of bacterial infections through the design of more specific and effective immune stimulants, adjuvants, and vaccines. Our results can also set the foundation for further developing of biomarkers, characterizing the mechanisms of immune organ barriers for vaccine development, and expanding our knowledge of teleost immunology.

Material and Methods

Ethics Statement

All experiments were performed according to local government regulations, and the procedures and protocols involving handling and treatment of fish were approved by the Institutional Animal Care and Use Committee of Qingdao Agricultural University.

Preparation of A. salmonicida subsp. masoucida Vaccines

Aeromonas salmonicida subsp. masoucida (strain RZ6S-1) was isolated from diseased cultured Atlantic salmon in a farm located in Yantai, Shandong province. The strain was stored in 15% glycerol (Sangon, China) at − 80 ℃ in laboratory until passage on tryptone soy agar supplemented with 1.5% NaCl (w/v) (TSA) for 50 h, immediately before use. Bacteria were cultured in tryptone soy broth supplemented with 1.5% NaCl (TSB) for approximately 12 h at 20 ℃ prior to use during the production of vaccine.

Aeromonas salmonicida subsp. masoucida (strain RZ6S-1) were cultured in TSB medium and centrifuged at 4000 g for 10 min. Bacteria were resuspended in 0.2 M phosphate buffer saline (PBS, 0.145 M Na2HPO4⋅12H2O, 0.055 M NaH2PO4⋅2H2O, pH 7.4), and inactivated with 0.2% formalin at 20 ℃ for 24 h, and 4 ℃ for 72 h. Success of inactivation of bacteria was checked by culturing on TSA for 7 days. Inactivated RZ6S-1 was suspended with PBS to approximately 2 × 109 CFU/mL and emulsified with the oil adjuvant Marcol 52 (MOBIL™, USA) at a ratio of 4:6 to form an oil-based vaccine (Yan et al. 2021).

Experimental Fish and Vaccine Injection

Atlantic salmon were obtained from a farm located in Yantai, Shandong province. Fish were cultured in ponds (7 m × 7 m × 0.5 m) equipped with recirculating aquaculture system (RAS) 14 ± 2 ℃, with continuous underground seawater (infiltrated coastal seawater, flow-through 6000–9000 L/h). The seawater salinity, the pH, and the dissolved oxygen was 2.7 ± 0.2%, 7.6 ± 0.1, and 7.0 ± 0.5 mg/L, respectively. Fish were fed twice (in case of disease outbreak only once) a day with commercial feed (Shengsuo, China), and the daily feeding rate was set at 1.2 ± 0.1% of fish weight at 0–8 weeks, 1 ± 0.1% at 9–16 weeks, and 0.9 ± 0.1% at 16–24 weeks. Before the vaccination, 20 fish were randomly selected, euthanized with an overdose of MS-222 (200 mg/L, Yibaolai, China), and examined microscopically for the presence of parasites. Samples were taken from internal organs for bacterial analysis on TSA plates to confirm that the fish were not infected by A. salmonicida.

A total of 5000 fish (33.56 ± 1.45 g) were randomly selected and equally divided into two ponds (2500 fish each), and then cultured with continuous underground seawater at 14 ± 2 ℃. One pond was designated as experimental group, and fish were intraperitoneally vaccinated using a repeater syringe with 0.1 mL Marcol 52 emulsified vaccines, containing 108 cells/fish of A. salmonicida subsp. masoucida. Fish of the control group were cultured in another pond under the same conditions without vaccination. All the fish were cultured as normal production by the experienced workers. The experiment was conducted over 3 months, and 30 fish (10 fish × 3 replicates) were randomly selected and euthanized with MS-222 (200 mg/L, Yibaolai, China) from each pond at 12 h (AS_12h), 24 h (AS_24h), 1 month (AS_1m), 2 months (AS_2m), and 3 months (AS_3 m) post-vaccination. The spleens from 10 fish were dissected and pooled together as one biological replicate, and a total of three replicate pooled samples were collected for each time point. All samples were flash-frozen in liquid nitrogen and then stored in an ultra-low freezer at − 80 ℃ until RNA extraction.

RNA Extraction, Library Construction, and Sequencing

Total RNA was extracted from the tissue samples using Trizol® Reagent (Invitrogen, USA), and 1% agarose gels were used to monitor the integrity of the RNA. After extraction, contaminating DNA was digested with RNase-Free DNase I (TIANGEN, Beijing, China) under the guidance of the instructions. Then, all RNA samples were run on the gel to ensure that no genomic DNA existed. The RNA quality and quantity of each sample were determined by Nanodrop 2000 (Thermo electronic North America LLC, FL). All extracted samples with A260/280 and A260/230 greater than 1.8 and RNA integrity number (RIN) greater than 7.0 were used for cDNA library construction. NanoPhotometer spectrophotometer (IMPLEN, CA, USA) was used to check RNA concentration. Three replicate samples were processed for each time point, and a total of 18 libraries were generated and sequenced.

NEBNext Ultra™ RNA Library Prep Kit for Illumina (NEB, USA) were used to generate sequencing libraries following manufacturer’s recommendations, and index codes were added to attribute sequences to each sample. Briefly, poly-T oligo-attached magnetic beads were used to purify mRNA from total RNA. Divalent cations were used to carry out fragmentation under elevated temperature in NEBNext® First Strand Synthesis Reaction Buffer (5X). First strand cDNA was synthesized using random hexamer primer and M-MuLV Reverse Transcriptase (RNaseH−). Subsequently, DNA polymerase I and RNase H were used to perform second strand cDNA synthesis. Remaining overhangs were converted into blunt ends via exonuclease/polymerase activities. After adenylation of 3′ ends of DNA fragments, NEBNext® Adaptors with hairpin loop structure were ligated to prepare for hybridization. The library fragments were purified with AMPure XP system (Beckman Coulter, Inc., Beverly, USA) to select cDNA fragments of preferentially 150–200 bp in length. Then, 3 µl USER Enzyme (NEB, USA) was used with size-selected, adaptor-ligated cDNA at 37 °C for 15 min followed by 5 min at 95 °C before PCR. Then, PCR was performed with Phusion® High-Fidelity DNA Polymerase, Universal PCR primers and Index (X) Primer. At last, PCR products were purified (AMPure XP system) and library quality was assessed on the Agilent 2100 Bioanalyzer system. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using TruSeq SR Cluster Kit v3-cBot-HS (Illumina) according to the manufacturer’s instructions. After cluster generation, eighteen spleen sequencing libraries (three replicates of 0 h, 12 h, 24 h, 1 month, 2 months, and 3 months, respectively) were sequenced on an Illumina Hiseq 2500 platform, and 150 bp paired-end reads were generated.

Reads Mapping and Differential Expression Analysis

FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) was used to assess the quality of sequencing data before and after trimming to make sure high-quality sequences were used in subsequent analysis. Trimmomatic (version 0.32) (Bolger et al. 2014) was used to trim raw reads by removing adapter sequences and ambiguous nucleotides. Only reads where both pairs had a length greater than 30 bp post-filtering were retained. Hisat2 was used to map the filtered reads to the Atlantic salmon genome reference (Lien et al. 2016), 5% mismatch of the mapped length was allowed and the minimum 90% of the bases matched to the genome was restricted. The HTSeq-count was used to extract read counts from the mapping profiles with the recommended mode (Anders et al. 2015).

EdgeR (Robinson et al. 2010) was used to conduct the trimmed mean of M-values normalization of expression value for each group in order to illustrate differentially expressed pattern. Differential expression between treatment and control samples was calculated with DEGseq (version 1.18.0) (Likun et al. 2010), and resulting P values were corrected for false discovery rate (FDR). DEGs were defined as showing FDR corrected P value < 0.05 and absolute value of log2 (fold change) > 1. Venn diagram and volcano plot was used to intuitively display the DEGs of different time points.

Gene Ontology and Enrichment Analysis

DEGs were annotated by gene ontology (GO) functional enrichment with the GO-TermFinder (V 0.86), and P values < 0.05 were considered significant difference (Boyle et al. 2004). The DEGs were enriched into three GO categories including biological process, molecular function, and cellular component at level 2. Pathway analysis was performed with the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg) (Kanehisa 2013), combined with manual literature review. Pathways with q-value < 0.05 were significantly enriched in DEGs after multiple test correction.

Experimental Validation by Quantitative Real-Time PCR

To confirm the reliability of RNA-seq results, we randomly selected 18 immune-related DEGs for qRT-PCR validation. The specific primers of all these genes were designed by PrimerQuest (https://sg.idtdna.com/PrimerQuest/Home) and listed in Table 1. The specificity of primers was assessed by aligning with the Atlantic salmon whole genome database (Lien et al. 2016) using BLASTN with E-value of 1e−10. Reverse transcription was performed with 1 μg total RNA using the qScript™ cDNA Synthesis Kit (Quanta Bioscience, Gaithersburg, MD). All the cDNA products were then diluted to 250 ng/μL with DNase/RNase-Free water (Solaribio, Beijing, China). The qPCR was performed using SYBR Green PCR Master Mix on a CFX96 real-time PCR detection system (Bio-Rad, USA). The cycling conditions of qPCR were denaturation, 95 °C/30 s, 40 cycles of 95 °C/5 s, 60 °C/5 s, and 72 °C/5 s. Test PCR was performed in advance to ensure all the genes were amplified with expect PCR product sizes. A no-template control was run on all plates. The Atlantic salmon β-actin gene was used as an internal reference gene (Ingerslev et al. 2006). For the transcriptome libraries preparation, first strand cDNA (500 ng RNA per 10 μL reaction) was synthesized by PrimeScript™ RT reagent Kit (Takara, Otsu, Japan) according to manufacturer’s protocol. The expression levels of those genes were detected by a CFX96 real-time PCR detection system (Bio-Rad Laboratories, Hercules, CA) following the manufacturer’s instructions. The reactions were performed in 20 μL volumes containing 2 μL of the diluted template cDNA, 0.6 μL of each primer (10 μM), 10 μL of SYBR® Premix Ex TaqTM II (TliRNaseH Plus), and 6.8 μL of RNA-free water. The thermal cycling profile was performed as follows. The PCR reaction mixture was denatured at 95 °C for 30 s and followed by 40 cycles of 95 °C for 30 s, 58 °C for 30 s and followed by dissociation curve analysis, 5 s at 65 °C, then up to 95 °C at a rate of 0.1 °C/s increment, to verify the specificity of the amplicons. No-template controls were performed as negative control for each gene, and each reaction was confirmed by repeating in triplicates for the qPCR analysis. Finally, 2−ΔΔCt method was used to calculate the relative gene expression fold changes (Livak and Schmittgen 2001). The qPCR was repeated in triplicate runs to confirm expression patterns.

Table 1

Primers used for qPCR validation

Gene symbol	Primers F	Primers R	Product size (bp)	Amplification efficiency	Accession no
Beta-actin	AAGATGAAATCGCCGCAC	ATGGAGGGGAAGACAGCC	61	108.61%	AF012125
IL8	CTCTTCTCATAGCGGCAGATTC	GTCAGTCCAGTTGCTAGTTCTT	32	110.15%	LOC100195681
IL1B-L	AGCGCAGTGGAAGAACATATAG	TCAGCACCCAGCACTTATTC	100	95.48%	LOC106570815
IL12B-L	AATGGTTCCTTCAACCGAGAG	CCTGTTCACCCATCATCCATAG	66	104.82%	LOC106612126
CCL11-L	TGTTACACTGGCGGAAGATG	GTGCTTGTCGTCCTTGTTTATTT	63	108.57%	LOC106564720
Pf2-L	GAGACAGAGAGCAGACGTATTG	GCTCAGAGTGGCAATGATCT	57	104.43%	LOC106601371
MIP2A	TAACACCATTAGACTTTATTAGAGG	GCCTCTCAGACTCATTCCC	92	100.47%	LOC100196393
H2-Aa-L	CCCATCAGAGACTGAAAGAGTTAG	AAACGGCTCCAGTGAGAATC	63	108.71%	LOC106567571
TNF-alpha-2	CTGGGCTCTTCTTCGTTTACA	GAGTCCGAATAGCGCCAAATA	64	108.74%	LOC100136458
TNFRST11B-L	GGCCAGTTACTCACCTGTAATC	CTGAGTGTAGTGGCTTGATGG	68	109.95%	LOC106596966
TNFRST10B-L	TGGACTGGCCTGAATGTATTT	TCCCAAATGCCACCCTATTT	79	102.09%	LOC106580524
PTGS2	CAAGATCACCTTTGTGCGAAAC	CCCAGCTCTTGTAGCCATAATC	66	99.86%	LOC100136456
FADD-L	CATTGGAAGGAAGTGGCAAAC	CAGGTTATTCGGGTGCTTCT	50	102.43%	LOC106561969
TNF-L	ACGGGCCTATACTTCGTCTAT	GTGTTACTCAGGTGAACCATCTC	57	106.84%	LOC106589058
PDGFC	CAGATCCGGGTGCGATTTAT	GAGTTACCTCGGGCTTCAATAC	64	107.69%	LOC106613208
CALCRL-L	GATGAACAGAGATCAGCCGTATAA	GGAGTCGAAGTCACGGAAATAG	93	95.58%	LOC106588943
CCL4-L	GGTACTCAGTGCCATCACTTT	TGGTGTAACCCACCACTTTAC	74	107.99%	LOC100136511
CCR9	CCTCAATCCCTTCCTGTATGTG	CCCTGGATCTTGCTGAGTTT	72	103.13%	LOC100195489
NFATC1-L	TCCTCCCGCCTATCCTTAAA	GTTGGGCTGACATATTGCTTTG	62	108.87%	LOC106569777

The complete name was shown in the supplementary file

Results

Sequencing and Mapping of Short Expressed Reads

RNA-seq analysis was performed to generate the expressed short reads from the spleen of Atlantic salmon in control group (0 h) and vaccination-treated groups (12 h, 24 h, 1 month, 2 months, and 3 months), and the results are summarized in Table 2. A total of 441,630,681 clean reads were obtained after removing low-quality reads, rRNA reads, reads containing adapters, and reads with > 10 unknown nucleotides. The sequence data were deposited at the NCBI Sequence Read Archive (SRA) with Bioproject number PRJNA786357. After filtering, the Q20 and Q30 were 94.86–96.87% and 88.44–92.57%, respectively, indicating that the sequencing data were of high quality. In addition, the average percentage of sequences successfully mapped onto the Atlantic salmon reference genome was 88.17%, resulted in a total of 389,373,457 mapped reads (Table 2).

Table 2

Summary of RNA-seq of Atlantic salmon spleen transcriptomes post-ASM vaccination

Sample name	Clean reads (bp)	Clean bases (bp)	Clean GC	Clean Q20	Clean Q30	Total mapped reads (%)
Spleen-0 h-1	23,990,624	6,517,245,648	48.26	96.58	92.04	21,406,834 (89.23%)
Spleen-0 h-2	21,784,367	7,232,369,024	48.15	96.73	92.34	19,468,689 (89.37%)
Spleen-0 h-3	24,176,459	8,241,674,628	48.57	96.87	92.57	21,282,537 (88.03%)
Spleen-12 h-1	27,647,896	7,577,562,270	48.31	95.69	89.97	24,108,966 (87.20%)
Spleen-12 h-2	25,369,710	8,589,785,570	48.08	94.88	88.45	21,894,060 (86.30%)
Spleen-12 h-3	28,744,767	7,158,238,314	47.45	94.86	88.45	25,263,776 (87.89%)
Spleen-24 h-1	23,979,715	6,212,895,586	48.58	94.93	88.56	20,975,057 (87.47%)
Spleen-24 h-2	20,939,499	7,208,112,598	46.44	95.68	89.97	18,805,765 (89.81%)
Spleen-24 h-3	24,232,981	6,028,183,378	47.31	95.05	88.82	21,475,268 (88.62%)
Spleen-1 m-1	20,190,911	6,392,119,708	47.98	96.23	91.34	17,832,613 (88.32%)
Spleen-1 m-2	21,396,364	6,386,269,348	46.93	94.89	88.44	18,963,598 (88.63%)
Spleen-1 m-3	21,428,207	6,943,521,486	47.11	95.81	90.50	18,878,251 (88.10%)
Spleen-2 m-1	23,197,668	6,802,255,556	48.57	96.32	91.49	20,516,018 (88.44%)
Spleen-2 m-2	22,837,050	8,049,544,156	46.32	95.47	89.60	20,299,854 (88.89%)
Spleen-2 m-3	26,922,069	9,479,971,644	48.30	96.79	92.55	23,594,502 (87.64%)
Spleen-3 m-1	31,850,655	8,593,044,178	47.98	96.86	92.44	28,334,343 (88.96%)
Spleen-3 m-2	28,887,121	7,133,708,518	47.21	95.88	90.72	25,677,762 (88.89%)
Spleen-3 m-3	24,054,618	6,517,245,648	48.47	95.57	89.62	20,595,564 (85.62%)
Total	441,630,681	131,063,747,258				389,373,457 (88.17%)

Analysis of Differentially Expressed Genes

In order to identify DEGs in Atlantic salmon spleen, pairwise comparisons were performed between control groups (0 h) and experimental groups (12 h, 24 h, 1 month, 2 months, and 3 months). The number of DEGs at different time points post-ASM vaccination are presented in Fig. 1 and Table 3. A total of 1872 DEGs were identified at 12 h post-vaccination, with 1125 DEGs up-regulated and 747 DEGs down-regulated. Subsequently, there were 4752 DEGs identified at 24 h post-vaccination, with 2126 DEGs up-regulated and 2626 DEGs down-regulated. Moreover, a total of 1916, 1074, and 2259 DEGs were identified at 1, 2, and 3 months post-vaccination, with 1098, 820, and 1351 DEGs up-regulated, and 818, 254, and 908 DEGs down-regulated, respectively (Fig. 1A–E and Table 3). Noteworthily, the number of up-regulated DEGs was higher than down-regulated DEGs at all tested time points, except for the 24 h. For the intersection analysis, a total of 119 DEGs were shared among all different time points post-vaccination (Fig. 1F).

Fig. 1

Analysis of DEGs displayed through volcano plot and Venn diagram in Atlantic salmon spleen at different time points post-vaccination. A–E Volcano plot of DEGs in Atlantic salmon spleen at 12 h, 24 h, 1 month, 2 months, and 3 months, respectively. The x-axis represents the value of log2 (fold change), and the y-axis represents the value of –log10 (FDR). Red and green spots represent significantly up-regulated and down-regulated genes, respectively (|log2 (fold change) |> 1) while blue spots represent no difference in gene expression. F Venn diagram

Table 3

Summary of DEGs in Atlantic salmon spleen at different time points post ASM vaccination with a cutoff of FDR-corrected P value < 0.05 and |log2 (fold change) |> 1

DEGs	12 h	24 h	1 month	2 months	3 months
Up-regulated	1125	2126	1098	820	1351
Down-regulated	747	2626	818	254	908
Total	1872	4752	1916	1074	2259

Analysis of Expression Signatures, Clustering, and GO Analysis

The expression signatures and hierarchical clustering of all DEGs were displayed with heatmap diagram. On the one hand, spleen samples of Atlantic salmon in control and vaccination groups at 2 months were firstly clustered together, and then clustered with vaccination-treated groups at 1 month and 3 months. On the other hand, the spleen samples of Atlantic salmon in 12 h and 24 h were clustered together, and finally fell into the other clade mentioned above. In addition, six major clusters of DEGs were observed based on their expression signatures (Fig. 2). For instance, the cluster 1 of DEGs were gradually up-regulated at different time points. Meanwhile, the cluster 6 of DEGs were firstly gradually down-regulated at 12 h and 24 h, and then gradually up-regulated at 1 month, 2 months, and 3 months.

Fig. 2

Heatmap display of hierarchical clustering of all DEGs in Atlantic salmon at different time points post-ASM vaccination. The DEGs are displayed in rows while samples are displayed in columns

To explore the biological functions of DEGs, GO analysis was performed at level two. Enriched GO terms were classified to biological process (BP), cellular component (CC), and molecular function (MF). On the one hand, the GO analysis results of up-regulated DEGs are presented at Fig. 3A. In BP, single-organism-metabolic-process (GO:0,044,723), immune system process (GO:0,002,376), and response to external stimulus (GO:0,009,605) were the most annotated terms. In MF category, oxidoreductase activity (GO:0,016,491) was the most annotated term. Under CC category, membrane (GO:0,005,886), membrane-part (GO:0,044,425), intrinsic component of membrane (GO:0,031,224), and integral component of membrane (GO:0,016,021) were the most common annotated terms. On the other hand, the GO analysis results of down-regulated DEGs are presented in Fig. 3B. In BP category, biological regulation (GO:0,065,007) and regulation of biological process (GO:0,050,789) were the two most enriched GO terms. Under molecular function category, signal transducer activity (GO:0,004,871), receptor activity (GO:0,004,872), and molecular transducer activity (GO:0,060,089) were the most enriched GO terms. Under cellular component category, extracellular region (GO:0,005,576) and extracellular matrix (GO:0,030,198) were the GO terms represented by the highest numbers of DEGs.

Fig. 3

Classification map of DEG secondary entries. The x-axis represents the second-level GO entry, and the y-axis represents the number of DEGs in the GO entry. A GO analysis of up-regulated DEGs. B GO analysis of down-regulated DEGs

Signaling Pathway Analysis

KEGG analysis significantly revealed a large number of immune-related pathways. A total of 2305 genes were assigned into KEGG functional pathways in Atlantic salmon spleen after vaccination treatment, and forty representative pathways were selected and presented in Fig. 4. On the one hand, 1050 up-regulated DEGs were assigned to twenty pathways such as cytokine-cytokine receptor interaction (116 DEGs), NOD-like receptor signaling pathway (40 DEGs), tumor necrosis factor (TNF) signaling pathway (61 DEGs), NF-kappa B signaling pathway (52 DEGs), intestinal immune network for lgA production (32 DEGs), and antigen processing and presentation (37 DEGs). On the other hand, 1255 down-regulated DEGs were enriched into another twenty pathways including pathways in cancer (169 DEGs), cytokine-cytokine receptor interaction (79 DEGs), intestinal immune network for lgA production (24 DEGs), PI3K-Ak signaling pathway (127 DEGs), and B cell receptor signaling pathway (41 DEGs).

Fig. 4

Pathway analysis of DEGs between control group and vaccination-treated groups based on KEGG database (FDR-corrected P value < 0.05). The vertical coordinates indicate the KEGG channel. The horizontal coordinates indicate the Rich factor. The larger the Rich factor, the greater the degree of enrichment. The larger the point, the greater the number of differential genes enriched in the pathway. A Pathway analysis of up-regulated DEGs. B Pathway analysis of down-regulated DEGs

Based on a combination of enrichment analysis, pathway analysis, manual annotation, and literature researches, representative DEGs were arranged into five categories according to functional pathways, including cytokine-cytokine receptor interaction, HTLV-I infection, MAPK signaling pathway, PI3K-Akt signaling pathway, and TNF signaling pathway (Table 4). Putative pathways involved in immunity response in Atlantic salmon spleen after vaccination treatment were illustrated in the context of a diagram (Fig. 5).

Table 4

Representative DEGs and their log2FC value in Atlantic salmon spleen post-ASM vaccination at different time points involved in key molecular pathways

Gene name	Log2 FC
	12 h	24 h	1 month	2 months	3 months
Cytokine-cytokine receptor interaction
C–C chemokine receptor type 7	−1.11	−1.15	−0.9	−0.93	−1.05
C–C chemokine receptor type 9	−0.09	−1.04	−1.39	−1.41	−2.33
C–C motif chemokine 4	4.5	4.67	2.84	3.05	4.03
CCL4-like chemokine	−0.59	−1.17	−1.1	−1.33	−1.41
C-X-C chemokine receptor type 3	−1.28	− 0.34	−1.6	−1.13	−2.57
C-X-C motif chemokine 11	2.17	1.47	1.85	1.51	1.23
Interleukin-17 receptor A	0.16	−0.76	−1.44	−1.17	−1.57
Tumor necrosis factor receptor superfamily member 10A	−0.64	−1.61	−1.4	−1.17	−1.27
Tumor necrosis factor receptor superfamily member 10B	1.02	1.63	1.43	0.44	1.36
Tumor necrosis factor receptor superfamily member 11B	7.65	7.95	5.8	5.81	6.51
Tumor necrosis factor receptor superfamily member 16	−1.13	−2.45	−1.77	−0.92	−1.12
Tumor necrosis factor receptor superfamily member 5	0.92	1.39	1.12	0.73	1.03
HTLV-1 infection
BOLA class I histocompatibility antigen, alpha chain BL3-7	1.26	1.67	1.87	1.81	2.03
Calcitonin gene-related peptide type 1 receptor	2.3	2.24	0.88	1.56	0.52
cAMP-responsive element modulator	5.57	4.8	2.05	2.03	1.96
Major histocompatibility complex class I-related gene protein	1.35	2.08	3.11	2.67	2.4
Platelet-derived growth factor receptor alpha-like	1.35	1.57	0.57	2.14	1.65
Proto-oncogene c-Fos	2.91	1.02	2.14	3.2	1.28
MAPK signaling pathway
Calcium voltage-gated channel auxiliary subunit alpha2delta 4	1.94	1.73	1.46	1.32	1.58
Interleukin-18 binding protein	−2.22	−2.56	−2.21	−0.89	−2.97
Interleukin-1 beta	6.53	3.9	3.3	3.06	1.85
Mitogen-activated protein kinase MLT	2.29	2.18	1.43	1.55	0.94
Tumor necrosis factor alpha-1 precursor	3.73	−0.82	4.42	5.18	4.32
Tumor necrosis factor alpha-2 precursor	7.66	4.33	7.35	8.23	6.71
PI3K-Akt signaling pathway
Collagen alpha-1(II) chain	8.14	-	-	4.36	-
H-2 class II histocompatibility antigen, A-K alpha chain	1.55	1.53	1.41	1.5	0.67
Phosphatidylinositol 3-kinase regulatory subunit gamma	1.17	1.4	1.62	2.25	1.43
Protein kinase N3	2.05	2.16	1.5	1.03	1.84
Thrombospondin-4-B	2.87	2.82	4.52	5.05	4.41
TNF signaling pathway
C-X-C motif chemokine 10	1	0.13	2	1.88	1.74
FAS-associated death domain protein	1.44	2.21	1.22	0.45	0.87
Matrix metalloproteinase-14	2.11	1.63	1.2	1.74	1.21
Prostaglandin G/H synthase 2	6.92	5.41	4.62	5.17	4.37
Proto-oncogene c-Fos	2.91	1.02	2.14	3.2	1.28
Transcription factor jun-B	4.82	3.83	3.1	2.58	2.62
Tumor necrosis factor receptor superfamily member 21	6.55	6.97	7.51	7.94	7.76

Fig. 5

Representative pathway analysis of DEGs between control group and vaccination-treated groups based on KEGG database (FDR-corrected P value < 0.05)

Validation of RNA-seq Results by qRT-PCR

To validate the DEGs identified by RNA-seq, we randomly selected eighteen immune-related genes for RT-qPCR validation. As shown in Fig. 6, fold changes at each time point from qPCR were compared with the RNA-seq expression profiles, and all trends of qPCR results were significantly correlated well with the RNA-seq expression profiles, suggesting that digital presentation of short reads in RNA-seq accurately reflects the transcriptome profiles.

Fig. 6

Validation of RNA-seq results using qPCR. The x-axis represents the time points, and the y-axis represents the relative fold change. The results of qPCR were presented as mean ± SE of fold change, and * indicates significance at the 0.05 level. The blue line indicated the variation trend of RNA-Seq, and the green square indicates significance at the 0.05 level

Discussion

Atlantic salmon is one of the most economically important cold-water fish in aquaculture in marine and freshwater ecosystems (Andrew et al. 2021; Wang et al. 2019), which is seriously harmed by A. salmonicida. A. salmonicida is one of the most important bacterial fish pathogens that threaten Atlantic salmon, causing huge losses to the global salmon farming industry (Dallaire-Dufresne et al. 2014; Janda and Abbott 2010). A. salmonicida have the capability to infect non‐phagocytic cells (Garduño et al. 2000; Valderrama et al. 2017) and phagocytic cells (Soto-Dávila et al. 2019). Once the fish is infected with the bacteria, the disease causes rapid septicemia resulting in the formation of necrotic lesions in the skin and hemorrhages in the internal organs, such infections are often fatal within as early as 2 or 3 days (Burr et al. 2005). Previous studies have revealed that A. salmonicida can greatly damage the immune organs of fish, after infection with A. salmonicida, symptoms such as exophthalmos, skin bleeding, ulcers, and necrosis in muscles and different internal organs (mainly spleen and kidneys) will occur (Burr et al. 2005). A. salmonicida subsp. masoucida (ASM) strains, recently isolated from several countries in Southeast Asia, were clustered in cluster VII among newly reported 14 clusters based on partial nucleotide sequences of the vapA gene (Gulla et al. 2016). Several ASM strains were reported to lead to huge economic damages to the local salmonid aquaculture in China (Du et al. 2015; Yan et al. 2021).

Vaccines are an effective way to prevent fish from being infected by bacterial diseases, and effective protection from A. salmonicida appeared to require injection of the vaccines (Krantz et al. 1963; Midtlyng et al. 1996b; Durbin et al. 1999). To contribute the knowledge base at the molecular level, we performed transcriptome analysis to identify the DEGs related to the Atlantic salmon spleen after intraperitoneal injection of the ASM vaccine, and better reflect the immune response of the spleen to the vaccination. In this study, a total of 441.63 million clean reads were obtained, and 389.37 million reads were mapped onto the Atlantic salmon reference genome. In addition, 1125, 2126, 1098, 820, and 1351 genes were significantly up-regulated, and 747, 2626, 818, 254, and 908 genes were significantly down-regulated post-ASM vaccination at 12 h, 24 h, 1 month, 2 months, and 3 months, respectively. It is worthy to mention that it would be more reasonable to set the time-matched mock-injection controls at each time point, since all comparisons in this study was conducted between pre-injected fish and immune-stimulated fish. The results of RT-qPCR showed a significant correlation with RNA-Seq results, suggesting the reliability of RNA-Seq for gene expression analysis. Many of the identified DEGs and overrepresented signaling pathways were involved in inflammation and immune response.

The innate immune system detects the presence and the nature of infection, provides the first line of host defense, and controls the initiation and determination of the effector class of the adaptive immune response (Zhang and Dong 2005). A specific immune response, such as the production of antibodies against a particular pathogen, is known as adaptive immune response, because it occurs during the lifetime of an individual as an adaptation to infection with that pathogen (Zhang and Dong 2005). Many signal transduction pathways participate in both innate and adaptive immune responses. From a combination of GO enrichment analysis and literature reviews, we can divide representative DEGs of Atlantic salmon into several major pathways such as cytokine-cytokine receptor interaction, HTLV-I infection, MAPK signaling pathway, PI3K-Akt signaling pathway, and TNF signaling pathway. Putative pathways involved in immune responses in Atlantic salmon spleen are illustrated as a diagram (Figs. 4 and 5). The putative functional roles and interaction of these signaling pathways are discussed below.

Cytokine-Cytokine Receptor Interaction

Cytokines are soluble extracellular proteins or glycoproteins that are crucial intercellular regulators and mobilizers of cells engaged in innate as well as adaptive inflammatory host defenses, cell growth, differentiation, cell death, angiogenesis, and development and repair processes aimed at the restoration of homeostasis. In this study, most of these DEGs enriched in this pathway were related to tumor necrosis factor receptors, interleukins, and chemokines. In detail, the TNF superfamilies are examples of signal transducers whose integrated actions impinge principally on the development, homeostasis, and adaptative responses of the immune system (Bodmer et al. 2002). In our work, the significant up-regulations of TNFRSF5, TNFRSF10b, and TNFRSF11b, and down-regulations of TNFRSF10a and TNFRSF16 were observed at almost all time points (Table 4). TNFRSF5 is constitutively expressed in a wide range of fish tissues, with relative high expression levels in immune tissues including spleen and head kidney (Park et al. 2005; Cai et al. 2017). It was reported that infection with Vibrio harveyi and i.p. injection of polyI:C could elevate TNFRSF5 transcript levels in the spleen of humphead snapper (Lutjanus sanguineus) (Cai et al. 2017). Noteworthily, DNA vaccination with viral hemorrhagic septicemia virus (VHSV) G-protein has been shown to up-regulate TNFRSF5 expression in Japanese flounder (Ju et al. 2005), which was consistent with our results. In addition, we found a number of chemokines and chemokine receptors differentially expressed after ASM vaccination at different time points. Chemokines promote leukocyte mobilization and regulate the immune responses and differentiation of the recruited cells (Alejo and Tafalla 2011; Esche et al. 2005). Besides their chemotactic function, several chemokines were reported to exert direct antimicrobial activity (Chan et al. 2008; Hasan et al. 2006; Nguyen and Vogel 2012). Increasing studies have been performing in aquaculture fish species, focused on immune function and their role in pathogenesis (Fu et al. 2017a, b, c; Wang et al. 2020; Zhao et al. 2020, 2021; Liao et al. 2018; Grimholt et al. 2015; Xu et al. 2020; Wang et al. 2020a, b; Li et al. 2021).

HTLV-1 Infection

Human T-lymphotropic virus type 1 (HTLV-1) is one of the most important members of the retrovirus family, and it infects peripheral-blood mononuclear cells (PBMC), particularly CD4 + T cells, leading to various complications due to changes in cytokines production and cell-signaling pathways (Futsch et al. 2018; Keikha et al. 2020). In the HTLV-1 infection pathway, MHC class I family, as a target for anti-tumor and allogeneic immunity, code for key proteins of the adaptive immune system, which present antigens from intra-cellular (MHC class I) and extra-cellular (MHC class II) pathogens (Minias et al. 2021; Ong et al. 2021). It plays a critical role in immunity by presenting peptides on the cell surface for T cell recognition (Mohan et al. 2021). In this study, MHC class I gene was significantly up-regulated with Log2FC of 1.35, 2.08, 3.11, 2.67, and 2.4 at 12 h, 24 h, 1 month, 2 months, and 3 months, respectively. In addition, platelet-derived growth factors (PDGFs) are among the best-characterized pro-fibrotic cytokines in systemic sclerosis, and it plays a key role in the pathogenesis of systemic sclerosis (Dragun et al. 2009). It has become apparent that members of the CREB family play important roles in the nuclear responses to a variety of external signals (Cesare et al. 1999). CREB and CREM play roles in many physiological systems, including memory and long-term potentiation, circadian rhythms, pituitary function, and spermatogenesis (Cesare et al. 1999). Through the sequencing results, it can be observed that the differential expression multiples of the CLEM gene in the HTLV pathway were the highest expression level at 12 h and 24 h post-vaccination, which indicates that this gene has a positive response to the vaccine. In this study, many immune-regulated genes in the HTLV-1 infection were identified as DEGs between vaccination-treated and control groups, suggesting a contribution of HTLV-1 infection to immune responses.

MAPK Signaling Pathway

The evolutionarily conserved mitogen-activated protein kinase (MAPK) cascades coordinately regulate a wide range of cellular processes, including cell proliferation, differentiation, metabolism, motility, survival, and apoptosis, and are therefore important for numerous physiological processes including innate and adaptive immune responses (Dickinson and Keyse 2006; Wagner and Nebreda 2009; Zhang and Dong 2005; Boon and Zhang 2016). Three major groups of MAPKs, namely the extracellular signal-regulated protein kinases (ERKs), the p38 MAP kinases, and the c-Jun NH2-terminal kinases (JNKs), have been identified in mammalian cells (Schaeffer and Weber 1999; Han and Ulevitch 1999; Davis 2000, Chang and Karin 2001). MAPK signaling responds to a broad range of extracellular signals from cytokines, growth factors, and environmental stresses (Krens et al. 2006; Tiedje et al. 2014). Meanwhile, MAPK signaling is subject to distinct spatiotemporal regulation by crosstalk and feedback mechanisms and can also cooperate with other pathways to control gene transcription (Menon and Gaestel 2016; Jordan et al. 2000). Therefore, the MAPK cascade acts as an important integrator of signal transduction, which is crucial for many physiological and pathological processes (Sun et al. 2015), such as cell proliferation, differentiation, apoptosis, cancer, and almost all aspects of immune responses (Kim and Choi 2015; Sun et al. 2015; Peti and Page 2013; Arthur and Ley 2013; Dong et al. 2002). Many DEGs related to immunity were found in the MAPK enrichment pathway, such as interleukin-1 beta (IL-1b), mitogen-activated protein kinase kinase kinase MLT-like (MAP3K20-L), tumor necrosis factor alpha-1 precursor (TNF-α-1), tumor necrosis factor alpha-2 precursor (TNF-α-2), calcium voltage-gated channel auxiliary subunit alpha2delta 4, and interleukin-18 binding protein. In detail, interleukin is a key proinflammatory cytokines that enables organisms to respond to infection and induces a cascade of reactions leading to inflammation (Dinarello 2002). MAP3K family is associated with aging, inflammation, oxidative stress, and related diseases, which plays an important role in animal innate immune response (Cao et al. 2019; Wang et al. 2020a, b). The cytokine TNF-α is a powerful proinflammatory cytokine released by several immune cells during infection or tissue damage and is involved in a diverse range of inflammatory, infectious, and malignant conditions (Roca et al. 2008). All of them indicate that the Atlantic salmon spleen might be an important organ in response to the vaccine. In addition, as a critical activator of the MAPK signaling pathway, map3k-l was induced in responses to vaccine with a log2FC of 0.94–2.29 in our study.

PI3K-Akt Signaling Pathway

The PI3K/Akt pathway is one of the major signaling pathways that have been identified as important players in regulating cell proliferation, growth, survival, and angiogenesis (Zhang et al. 2021). Activated PI3K/AKT pathway contributes to cancer cell proliferation, survival, angiogenesis and inhibition of apoptosis, which are important factors of tumorigenesis (Chen et al. 2017). Apart from regulation of tumors, current evidence suggests PI3K/AKT can also modulate inflammation, apoptosis and immune cells proliferation, differentiation through regulating the activity of downstream effector molecules, such as NF-κB, Bad, caspase-9, and m-TOR (Vivanco and Sawyers 2002; Way et al. 2016; Park et al. 2014). Additionally, numerous studies have indicated that the PI3K/AKT signaling pathway is an important intracellular signaling pathway that plays a role in a variety of cellular physiological processes, including cell proliferation, cell apoptosis, and inflammatory response (Li et al. 2013; Jiang et al. 2019; Wang et al. 2020a, b). Chang et al. reported that 7‑ketocholesterol contributed to thrombosis via the induction of endothelial damage, apoptosis, and inflammatory responses, which were associated with the activation of PI3K/AKT signaling (Chang et al. 2016).

In this study, thrombospondin-4-B-like (thbs4b-l), phosphatidylinositol 3-kinase regulatory subunit gamma (PIK3R3), protein kinase N3 (PKN3), collagen alpha-1(II) chain-like (Col2a1-l), H-2 class II histocompatibility antigen (H-2), A-K alpha chain-like genes related to this pathway were enriched. In detail, thrombospondin-4 (THBS4) is a member of the extracellular calcium-binding protein family. It is a secreted pentameric globular protein that forms part of the extracellular matrix, and functions in calcium binding, cell attachment, cell migration, cellular proliferation, cytoskeletal organization, neurite growth, binding of other extracellular matrix components, and cell to cell interactions (Stenina et al. 2003; Adams and Lawler 2004; Arber and Caroni 1995; Carlson et al. 2008; Si et al. 2003; van Doorn et al. 2005). The PIK3R3 regulatory subunit is crucial in cell proliferation, cell differentiation, tumor angiogenesis, and tumor growth (Wang et al. 2013, 2014, 2015). Protein kinase N (PKN) is a serine/threonine protein kinase with a catalytic domain homologous to protein kinase C and a unique regulatory region containing antiparallel coiled-coil (ACC) domains (Maesaki et al. 1999). PKN is composed of three isoforms (PKN1, PKN2, and PKN3) derived from different genes in mammals (Mukai et al. 2016). Previous reports have shown that ablation of the PKN3 gene suppressed migration of embryonic fibroblast cells induced by various growth factors (Mukai et al. 2016). So far, there have been accumulated reports about the potential function of PKN isoforms using cultured cell experiments such as involvement in the regulation of cytoskeletal reorganization (Vincent and Settleman 1997; Mukai et al. 1997), cell adhesion (Calautti et al. 2002, Wallace et al. 2011), cell-cycle regulation (Isagawa et al. 2005; Schmidt et al. 2007), and tumorigenesis (Metzger et al. 2003, Leenders et al. 2004). COL2A1 encodes for type II collagen, the most important and abundant extracellular matrix protein in epiphysis cartilage (Zhou et al. 2021). It is closely related to the development and maturation of chondrocytes, synthesis and catabolism of extracellular matrix, and cartilage ossification (Barat-Houari et al. 2016). Interestingly, the differential expression of COL2A1-1 is particularly up-regulated at 12 h after the injection of the vaccine, while there is almost no differential expression at 24 h and 1 month and 3 months after the injection. Its function in the spleen needs to be further studied.

TNF Signaling Pathway

TNF is a member of a superfamily of cytokines that consists of more than 20 structurally related transmembrane and soluble proteins, which play a crucial role in various biological events by recruiting several intracellular adaptors, thus activating multiple signal transduction pathways (Wajant et al. 2003). The proteins of the TNF receptor superfamily are a group of cell-surface receptors critically involved in the maintenance of homeostasis of the immune system (Tong et al. 2002). The members of the TNF receptor superfamily are intimately involved in the regulation of the proliferation and death of immune cells and are of particular interest in relation to their role in immuno-biology and the pathogenesis of autoimmune disease (Chan et al. 2000). Functionally, the members of the TNF receptor superfamily can deliver signals leading to apoptosis or proliferation/survival (Tong et al. 2002). In the present study, many DEGs related to tumor necrosis factor are enriched in this pathway, such as prostaglandin G/H synthase 2 (PTGS2), tumor necrosis factor receptor superfamily member 21-like (TNFRSF21-l), and C-X-C motif chemokine 10 (CXCL10). PTGS2 promotes inflammation and suppresses T-cell-mediated adaptive immunity (Cao et al. 2016). CXCL10 may play a role in chronic inflammatory diseases, including coronary artery disease and related manifestations of atherosclerosis (Heller et al. 2006).

Below is the link to the electronic supplementary material.

Supplementary file1 (DOCX 18 KB)