Investigation of the Transcriptomic Response of Atlantic Salmon (Salmo Salar) Exposed to Paramoeba Perurans Using Next Generation Sequencing

Amoebic Gill Disease (AGD), caused by the protozoan extracellular parasite Paramoeba perurans, is a disease affecting Atlantic salmon (Salmo salar) aquaculture. Many studies to date have investigated the pathogenesis of ADG focusing on the host immune response in the gill after the appearance of clinical symptoms. This study investigated the gill transcriptomic prole of pre-clinical AGD using RNA-sequencing (RNA-seq) technology. RNA-seq libraries generated at 4, 7, 14 and 16 days post inoculation (dpi) identied 29,357 Differentially Expressed Genes (DEGs). RNA-seq data was validated using real-time, quantitative PCR (qPCR) analysis of 10 selected immune genes. DEGs mapped to 224 Gene Ontology (GO) terms, 27 reference pathways in the Kyoto Encyclopedia of Genes and Genomes (KEGG) and 15 Reactome Gene Sets. Immune suppression was evident at 7 dpi, prior to there being any evidence of ADG on the gill, involving signalling pathways for interleukins, Nod-like receptors, B-cell and T-cell receptors, and the differentiation of Th1/Th2 and Th17 cells. The results of this study suggest a mechanism for how N. perurans circumvents the host immune response to establish a successful infection, and could potentially lead to the development of novel strategies for AGD mitigation or prevention in aquaculture.


Introduction
Amoebic Gill Disease (AGD) is a proliferative gill disease of marine cultured Atlantic salmon, with the aetiological agent being the free-living protozoan Paramoeba perurans 1 . Among the parasitic conditions affecting gill health, AGD is the most signi cant in terms of prevalence and economic impact 2 . First described in Tasmania in 1985 1 , AGD is now present in all Atlantic salmon producing countries. Infection of the gill with Paramoeba sp. is reported to induce a tri-phasic host response, including a localised reaction to adhered parasites, a non-speci c immuno-regulatory cell in ltration and advanced hyperplasia with epithelial strati cation laments 3 . Macroscopically, AGD lesions are visible as white raised mucoid patches on the gill, as a result of increased mucus production by mucous cells, which forms the basis of the gross gill score scheme developed by Taylor et al. (2009) 4 to monitor the progress of AGD. In advance of gill pathology, the presence of P. perurans can be detected using a diagnostic real-time PCR (qPCR) 5 . The treatment for diagnosed AGD is either freshwater (< 3 PSU for 3 hours) or hydrogen peroxide (H 2 O 2 ) (800-1300 ppm for < 20 minutes) baths 6 .
Previous AGD studies identi ed up-regulation of the pro-in ammatory cytokine, interleukin-1beta (IL-1b) as being the hallmark of advanced stage AGD in Atlantic salmon 11,13 . A recent study using 2D quantitative RT-PCR found AGD-affected gills displayed an increased mRNA expression of cellular markers of immune cells, including professional antigen presenting cells (MHC-II, CD4), B cells (IgM, IgT, MHC-II) and T cells (TCR, CD4, CD8) 14 .
The down-regulation of tumour suppressor p53 (p53) has been suggested as one of the mechanisms underlying cell proliferation in AGD 7 . The differential expression of genes involved in oxidative stress has also been reported 17 .
Th1, Th17 and T-reg pathways were found to be down-regulated in AGD infections initiated with a higher dose of amoeba (5000 amoeba/L) with the Th2 pathway up-regulated by both high and low infection doses (500 and 5000 amoeba /L) 10 . Th2 cytokines, (il4/13a and il4/13b2), and the mucin Muc5AC have also been found to be up-regulated in the gill in AGD 12 . The aim of the current study was to investigate the early-phase transcriptomic host response of naïve salmon inoculated with P. perurans using Next Generation Sequencing (NGS), speci cally RNA-seq, prior to the onset of gill pathology. As RNA-seq can be used as a discovery-based technique, it does not require prior knowledge of the genomic information related to genes of interest, with the advantage of facilitating the identi cation of both known and unknown transcripts.

Results
Clinical symptoms of AGD were determined by macroscopic examination of the intact gills and were scored ranging from 0 (absence of clinical symptoms) to 5 (extensive lesions covering most of the gills surface) 4 . No lesions were evident on the gills in AGD-infected groups at 4 and 7 dpi, with 3/6 sh at gill score 1 at 14 dpi, and 6/6 sh with gill score 1 at 16 dpi. Diagnostic qPCR 5 analysis of the gill arch samples con rmed the presence of P. perurans in 5/6 at 4 dpi, 6/6 sh at 7 dpi and 14 dpi and in 5/6 sh at 16 dpi. There was no P. perurans detected in any of the control sh sampled at any time point. Gills were examined microscopically following staining with haematoxylin/eosin at all time points (Fig. 2).
Evidence of hyperplasia and fusion of the lamellar epithelium ( Fig. 2b) was observed at 16 dpi.
A total of 29,357 genes were found to be differentially regulated in the gill as a result of AGD infection, with approximately equal numbers of up-and down-regulated genes ( Table 1). The %DEGs at 4, 14 and 16 dpi ranged from 16.5-21.9% with a higher % of genes down-regulated genes than up-regulated (61.7%, 52.4% and 54.7%, respectively). At 7 dpi there were 42.3% DEGs, with more up-regulated than downregulated (59.4% vs 40.4%).  Table 2).
At 14 dpi, the reactome gene set R-DRE-983695: antigen activates b-cell receptor (BCR) leading to generation of second messengers pathway was the only enriched pathway (  At 4, 14 and 16 dpi process and pathway enrichment identi ed 14 up-regulated genes with involvement in immune regulation and activation (Table 3). Immune biological processes and a reactome gene set ( Table 2) were identi ed from the DEG lists generated at each time point and are listed in Table 3 and included GO:0002682: Regulation of immune system process (btk, c1cq, CD79b, CD99, cgas, hmgb1a, hmgb3b, Il-34, kita, lyn, tfpi1, themis), GO:0002253: Activation of immune system process (btk, c1cq, CD79b, cgas, hmgb1a, themis), GO:0097190: B-cell differentiation (kita, cd79b and Ikzf1.7), GO:0002768: immune response-regulating cell surface receptor signalling pathway (kita, btk, cd79b, themis2). The reactome gene set is R-DRE-983695: Antigen activates B Cell Receptor (BCR) leading to generation of second messengers (dapp, btk, and CD79b). At 7 dpi the reactome gene set R-DRE-168256: immune system ( Table 2) was associated with 51 downregulated genes (Table 4). Within this reactome gene set, 16 genes were associated with R-DRE-1280215: cytokine signalling in the immune system with all but 2 of these genes associated with R-DRE-449147: signalling by interleukins ( Table 2). Screening of the DEG lists at 7 dpi identi ed multiple transcripts for 9 members of the interleukin family as being down-regulated (IL-1β, IL-8, IL-11, IL-12β, IL-15, IL-17D, IL-17F, IL-18, IL-34) with only one member, IL-27β signi cantly up-regulated (12.4 log 2 FC) (supplementary Table   S3 online). Functional enrichment analysis (string-de.org) using the 51 genes identi ed in R-DRE-168256: immune system genes at 7 dpi ( The pattern and extent of DEGs at 7 dpi was different from the other time points with more genes showing differential expression (42%) and with more genes being up-regulated (59.6%) ( Table 1).
Analysis of the DEG lists at 7 dpi found that of the top 100 up-regulated genes there were 87 characterised and 69 unique transcript (supplementary Table S4  The top 100 down-regulated genes at 7 dpi yielded 78 unique transcripts and 19 uncharacterised genes (supplementary Table S5 online). The genes that were represented by multiple transcripts included C-C motif chemokine 4-like (2x), collagenase 3-like (2x), interferon-induced very large GTPase 1-like (2x), NADPH oxidase organizer 1-like (2x). The gene with the most down-regulated expression was interferoninduced guanylate-binding protein 1-like (log 2 FC, -7.5). Multiple members of the same gene family included the c-c chemokine family (cc4, cc20), the interleukins (il-1b, Il-17F, il-17 receptor E), and mucins (muc2, muc7).   Kita is a tyrosine-protein kinase that acts as cell-surface receptor for stem cell factor (Scf) and plays an essential role in the regulation of cell survival and proliferation, hematopoiesis, stem cell maintenance, gametogenesis, mast cell development, migration and function, and in melanogenesis 19 .
Themis2 is a gene which encodes a protein that plays a regulatory role in both positive and negative Tcell selection during late thymocyte development. The protein functions through T-cell antigen receptor (TCR) signalling, and is necessary for proper lineage commitment and maturation of T-cells. Themis2 plays a role in the in macrophage in ammatory response, promoting LPS-induced TNF production.
Five genes were found to be up-regulated only at 14 dpi: CD99, high mobility group b3b (hmgb3b) interleukin-34 (IL-34), Lyn and tissue factor pathway inhibitor 2 (t p2). CD99 has been described as a costimulatory molecule on T cells 20  MAPKs are a superfamily of serine/threonine protein kinases that transduce a variety of external signals, leading to an array of cellular responses that include growth, differentiation, apoptosis, and host defence response.
The up-regulation of Lyn is only evident from 14 dpi, when mucoid patches were rst identi ed on the gills. Tfpi2 is a serine protease inhibitor and is thought to play a role in the regulation of plasminmediated matrix remodelling 26  C1qC was up-regulated suggesting the classical complement pathway was activated. Composed of C1, C4, and C2 components reacting in this order, the classical pathway primarily recognizes antibodies in immune complexes. However, an increased expression of immunoglobulins were not observed in this study indicating that C1q is potentially interacting with other acute phase molecules.
Cgas is a cytosolic DNA sensor that activates a type-I interferon response 29 , DAPP is a B-cell-associated adapter that regulates B-cell antigen receptor (BCR) 30  IL-27β, was the only interleukin signi cantly induced at 7 dpi. IL-27 is composed of two non-covalently linked subunits, IL-27p28 (p28) and IL-27β, also called Epstein-Barr-virus-induced molecule 3 (EBI3) 36 .
These subunits exhibit structural and sequence homology to IL-12 subunits and IL-6 37 . IL-27 is unique in that although it induces Th1 differentiation, the same cytokine suppresses immune responses. IL-27 can antagonise the development of the Th17-cell response and limit Th-17 driven in ammation 38 which are critical for host defence against bacterial, fungal and viral infections at mucosal surfaces 39 .
Of the down-regulated immune genes in reactome R-DRE-168256 at 7 dpi (Table 4), 8 genes were found to participate in the Nod-like receptor signalling pathway. Nod-like receptors (NLRs) can initiate or regulate host defence pathways through formation of signalling platforms that subsequently trigger the activation of in ammatory caspases and NF-kB 40 . Genes found to participate in the Nod-like receptor signalling pathway included cxcl8, hsp90AB1, ikbkb, irf9, mapk3, nlrx1, sugt1, and tbk1. Nod-like receptors (NLRs) sense pathogen-associated molecular patterns (PAMPs) (pathogens/foreign) or damage-associated molecular patterns (cells/self) 41 . lrx1 is a dsRNA receptor 42 and was identi ed as being down-regulated in our data set at all 4 time points (4, 7, 14 and 16 dpi).
Interestingly, all of these signalling pathways had 2 genes in common: inhibitor of nuclear factor kappa B kinase subunit beta (ikbkb) and mitogen-activated protein kinase 3 (mapk3).
The transcription factor Nuclear Factor-kappa beta (NF-κβ), prior to activation, is held in the cytoplasm by the attachment of inhibitor kappa beta (Iκβ) and the formation of the Iκβ/NF-κβ complex Ikbkb is a gene which encodes the enzyme Iκβ kinase beta (IKKβ) which can phosphorylate 2 serine residues on the inhibitor in the Iκβ/NF-κβ complex, leading to the dissociation and degradation of the Iκβ inhibitor and the subsequent activation of NF-κβ 43 . The dissociated NF-κβ can translocate into the nucleus and activate the transcription of hundreds of genes involved in immune response, growth control, or protection against apoptosis. IKKβ is critical for cytokine production via NF-kB activation. The initial innate immune response is under the control of IKKβ and culminates in a successful humoral response that is dependent on IKKα 44 . Mapk3 acts upstream of IKKβ in the canonical NF-κβ activation pathway 45 .
This is the rst study exploring the transcriptomic response of Atlantic salmon to P. perurans in a controlled environment. The data presented show that the host response of Atlantic salmon is activated in advance of any clinical symptoms developing on gill tissue of sh inoculated with P. perurans. Of particular interest is the immune suppression brought about through the downregulation of 2 key genes, ikbkb and mapk3/ERK1 resulting in the continued inhibition of NF-κβ by Iκβ in the cytoplasm. Pathogens have previously been reported to developed strategies to circumvent the activation of the NF-κβ activation, by preventing the inhibitor, Iκβ, from being ubiquitinated and therefore preventing its degradation, causing NF-κβ to remain sequestered in the cell cytoplasm and therefore inactive 46 . Indeed, some viruses encode virulence factors, for example vaccinia viral protein B14, that target IKKβ to inhibit NF-κβ-mediated antiviral immune response 47 , suggesting that virulence factors associated with P. perurans potentially have some of immunomodulatory effect on their host. The present study provides the initial discovery and description of genes showing differential expression during early-phase ADG exposure/infection, which provides the basis for future, more in-depth studies of AGD-related immune response pathways in Atlantic salmon.

Fish husbandry
Post Paramoeba peruransisolation and culture P. perurans trophozoites were collected by gill swabbing from AGD infected Atlantic salmon on a commercial Irish farm. Amoebae were cultured on marine yeast agar plates (MYA; 0.01% malt, 0.01% yeast, 2% Bacto Agar), 16°C overlaid with 7 ml sterile seawater 49 , and sub-cultured weekly by transferring free-oating cells to fresh MYA plates. Con rmation of P. perurans identity was performed using qPCR as previously described by Downes (2015) 5 .

Paramoeba peruranschallenge
After an acclimatisation period of 6 days, 90 sh (45 x 2) were challenged with P. perurans (2750 amoebae/L) in 300 L for 4 h with oxygen saturation, sh behaviour and welfare closely monitored. Controls, 90 sh (45 x 2) were also held at 300 L for 4 h. Following holding in challenge or control tanks, sh were placed into new tanks containing 1000 L seawater.

Disease progression
Disease progression was monitored using gill scoring carried out on euthanised sh (Taylor, 2009). Whole gill samples were taken for qPCR 5 and histology to con rm the presence of P. perurans. For histological analyses, samples were xed in 10% neutral buffered formalin, routinely processed, embedded in para n wax blocks, sectioned (3-5 µm), and stained with haematoxylin and eosin (H&E).
Total RNA extraction Differential expression analysis Differential expression analysis was performed on biological replicates (n = 6) from 4, 7, 14 and 16 dpi using the DESeq2 (Version 1.24.0) where each time point was compared to time 0, pre-AGD samples ( Fig. 1). DESeq R provide statistical routines for determining differential expression in digital gene expression data using a model based on the negative binomial distribution. The resulting P-values were adjusted using the Benjamin and Hochberg's approach for controlling the false discovery rate. Genes with an adjusted p-value (FDR adjusted) < 0.05 found by DESeq R were assigned as differentially expressed. Eight gene lists were generated: up-and down-regulated genes at 4, 7, 14 and 16 dpi.
Protein-Protein Interaction (PPI) PPI was carried out using STRING (version 11.0, https://string-db.org) on the 8 gene lists generated from the pairwise comparisons (4, 7, 14 and 16 dpi, up-and down-regulated genes) using the default settings and Danio rerio selected as the species of interest. Signi cant pathway enrichments and functional information for genes identi ed in the Immune System (R-DRE-168256) was carried out using STRING against the Homo sapiens database.

Pathway and Process Enrichment
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway mapping was performed using Metascape (metascape.org). The zebra sh (Danio rerio) database was used to determine GO enrichment. Pathway and process enrichment analysis was performed using the following ontology sources: KEGG Functional Sets, KEGG Pathway, GO Biological Processes, GO Cellular Components, GO Molecular Functions, KEGG Structural Complexes and Reactome Gene Sets.
A user-supplied list of 5824 genes was used as the enrichment background. Terms with a p-value < 0.01, a minimum count of 3, and an enrichment factor > 1.5 are collected and grouped into clusters based on their membership similarities.
Reverse transcription was carried out using the GoScript (Promega) kit as per manufacturer's instructions. A pre-ampli cation step was adopted for the multiplex ampli cation of the target genes using a MiniAmp Plus PCR machine (Applied Biosystems), as per manufacturer's recommended protocols (PN 100-5875 C1, Fluidigm). The pre-ampli ed cDNA was treated with Exonuclease I to remove unincorporated primers prior to running on a Biomark HD, as per manufacturer's recommended protocols (PN 100-9791 B1, Fluidigm). PCR ampli cation e ciency (E) was calculated for each gene of interest and the housekeeping gene by the generation of standard curves using 10-fold serial dilutions of the cDNA template (standard curve: R 2 > 0.980, ampli cation e ciency range 90-105%). Melt curve analysis was used to con rm the ampli cation of single, PCR products.

Statistical analysis
For pathway and process enrichment analyses, p-values < 0.01 were calculated based on the accumulative hypergeometric distribution, and q-values calculated using the Benjamin-Hochberg correction for multiple testing using Metascape (metascape.org). Kappa scores were used as the similarity metric when performing hierarchical clustering on the enriched terms, and sub-trees with a similarity of > 0.3 were considered a cluster. The most statistically signi cant term within a cluster was chosen to represent the cluster. For qPCR analysis, the fold change of each gene at each time point was analysed relative to the T0 control using an un-paired t-test with differences considered signi cant at p < 0.05. Ggplot2_3.2.1 in R studio Version 3.5.1 was used to generate the Spearman correlation data. Plotly.py version 4.0.0 in Python 3.7.3 was used to graph the RNA-seq and qPCR data validation data.

Declarations Data Availability Statement
All data supporting this study are included in the results section, the supplementary material section or openly available in public databases.