Transcriptomic Profiling in Fins of Atlantic Salmon Parasitized with Sea Lice: Evidence for an Early Imbalance Between Chalimus-Induced Immunomodulation and the Host's Defense Response.

Parasitic sea lice (e.g., Lepeophtheirus salmonis) cause costly outbreaks in salmon farming. Molecular insights into parasite-induced host responses will provide the basis for improved management strategies. We investigated the early transcriptomic responses in pelvic fins of Atlantic salmon parasitized with chalimus I stage sea lice. Fin samples collected from non-infected (i.e. pre-infected) control (PRE) and at chalimus-attachment sites (ATT) and adjacent to chalimus-attachment sites (ADJ) from infected fish were used in profiling global gene expression using 44 K microarrays. We identified 6568 differentially expressed probes (DEPs, FDR < 5%) that included 1928 shared DEPs between ATT and ADJ compared to PRE. The ATT versus ADJ comparison revealed 90 DEPs, all of which were upregulated in ATT samples. Gene ontology/pathway term network analyses revealed profound changes in physiological processes, including extracellular matrix (ECM) degradation, tissue repair/remodeling and wound healing, immunity and defense, chemotaxis and signaling, antiviral response, and redox homeostasis in infected fins. The QPCR analysis of 37 microarray-identified transcripts representing these functional themes served to confirm the microarray results with a significant positive correlation (p < 0.0001). Most immune/defense-relevant transcripts were downregulated in both ATT and ADJ sites compared to PRE, suggesting that chalimus exerts immunosuppressive effects in the salmon's fins. The comparison between ATT and ADJ sites demonstrated the upregulation of a suite of immune-relevant transcripts, evidencing the salmon's attempt to mount an anti-lice response. We hypothesize that an imbalance between immunomodulation caused by chalimus during the early phase of infection and weak defense response manifested by Atlantic salmon makes it a susceptible host for L. salmonis.

salmon, respectively; whereas, no DEGs were discovered in Atlantic salmon fins [8]. Consequently, the mechanisms related to local responses in fins, which are the initial and preferential settlement sites of chalimi, against the early phase of L. salmonis infection are not completely understood. The present study addressed this gap of knowledge in the salmon's anti-parasitic immune response using a 44K consortium for Genomic Research on All Salmonids Project (cGRASP) microarray platform [27] which provides good coverage of the salmon transcriptome, together with real-time quantitative polymerase chain reaction (QPCR) for confirmation purposes. We have used this approach in our previous immune-relevant experiments for studying salmon immune transcriptome responses against viral mimic [28] and bacterial antigens [29]. We investigated the initial local host transcriptomic response in pelvic fin, the first dominant attachment sites of L. salmonis copepodids [30], of infected salmon by profiling the gene expression at attachment sites (ATT) and adjacent to chalimus-attachment sites (ADJ) compared to non-infected controls (PRE). Our study was aimed at depicting the modulation of host responses against L. salmonis infection, which enabled us to explore, in detail, specific genes, pathways, and networks associated with the interface of host-parasite interaction.

Sea Lice Load in Salmon
In this study, we used a host-parasite (S. salar-L. salmonis) model to investigate the transcriptomic responses of the host to infection. Pelvic fin tissue was sampled from non-infected fish (PRE; 0 days post-infection (dpi); n = 11), and salmon were experimentally infected with L. salmonis copepodids. At 8 dpi, the number of chalimi on each infected salmon was counted prior to the excision of fin tips for sampling from chalimus-attachment (ATT) and adjacent (ADJ) sites of each fish (n = 12).

Global Transcriptomic Changes in Response to Chalimus Attachment
In order to explore the transcriptional response in fin induced by the attached L. salmonis chalimi on pelvic fins during the infection (at 8 dpi), the differentially expressed probes (DEPs) between fins from PRE, ADJ, and ATT groups (n = 6 samples/group) were identified using a 44K cGRASP salmon microarray platform [27].
The fin transcriptome was compared between sample groups using Significance Analysis of Microarrays (SAM) [31], and the profile of DEPs has been illustrated in Figure 1A,B. The current study identified a total of 6568 non-redundant DEPs when three different groups were compared in pairs (i.e., 6227 in ADJ compared to PRE (ADJ versus PRE), 2239 in ATT compared to PRE (ATT versus PRE), and 90 in ATT compared to ADJ (ATT versus ADJ)) at false discovery rate (FDR) < 5% ( Figure 1A; Supplementary File S1). Of these DEPs, putative identities could be determined for 5373 DEPs by using BLASTn/x searches against NCBI nr/nt databases. A more stringent FDR cutoff of 1% yielded 1377, 677, and 42 chalimus-responsive transcripts for ADJ versus PRE, ATT versus PRE, and ATT versus ADJ comparisons, respectively ( Figure 1A). All downstream analyses presented here are based on transcripts identified at FDR < 5%.
Among 6227 DEPs in the ADJ versus PRE comparison, 3817 probes were significantly upregulated and 2410 probes were significantly downregulated ( Figure 1B). We identified 779 higher and 1460 lower expressed probes in the ATT versus PRE comparison, respectively. All 90 DEPs in the ATT versus ADJ comparison were significantly upregulated in ATT. The re-annotation of these probe lists provided information on DEGs for each comparison: 4260 DEGs for ADJ versus PRE, 1578 DEGs for ATT versus PRE, and 70 DEGs for ATT versus ADJ ( Figure 1B). By taking the redundancy of the 44K platform into account, it was determined that 6568 DEPs represent 3571 DEGs with known putative identity. For complete information of these DEPs, including their identity based on BLAST annotation and fold-change (FC) values, refer to Supplementary Files S1 and S2.
A relatively clear clustering of PRE, ATT, and ADJ samples was evident from the principal coordinate analyses (PCoA) and hierarchical clustering analyses conducted based on gene expression data of the entire set of DEPs (6568; Figure 1C,D). PCO1 and PCO2 accounted for 62.1% and 10.0% of the variation among fish, respectively. Along with PCO1, animals were distinguished according to the treatment, with samples obtained from the PRE control clustered distantly from samples obtained from infected fish (ATT and ADJ; Figure 1C). Moreover, the grouping pattern in hierarchical clustering ( Figure 1D) was also in good agreement with the PCoA results ( Figure 1C), in which ATT and ADJ samples clustered closely with considerable overlap. These results suggest that ATT and ADJ samples (both collected from the same fish) may share an overall gene expression pattern.   (1)(2)(3)(4). Bold numbers indicate the DEPs contributing to a particular segment. The number of DEGs corresponding to DEPs of each segment is shown. FC_# indicates the number of DEPs with a fold-change ≥ |2|. It should be noted that DEPs in segment 1 required to comply with the criteria of FC ≥ |2| in both comparisons ADJ versus PRE and ATT versus PRE. Only the features with FC ≥ |2| were subjected to Blast2GO enrichment analyses, except for segment 2 as this probe list was relatively shorter. Whereas, all the DEGs of each segment were used in gene ontology (GO)/pathway terms analyses using ClueGO. vs., versus.

Enriched Gene Ontology (GO) Terms Associated with DEPs
By increasing the FC cutoff to |2|, the DEPs were filtered to identify those likely to be biologically most relevant in terms of response to sea lice. This reduced the number of chalimus-responsive probes to 886 DEPs for ADJ versus PRE (219 upregulated and 667 downregulated), 537 for ATT versus PRE (108 upregulated and 429 downregulated), and 57 for ATT versus ADJ (all upregulated) ( Figure 1B, Supplementary File S1). Gene ontology (GO) term analyses of DEPs were conducted by Blast2GO in two steps. In the first step, GO term distribution of different GO domains (biological process (BP), molecular function (MF), cellular component (CC)) was mapped. Secondly, Fisher's Exact Test was conducted to study the enrichment of GO terms. Figure 3 illustrates the results of these two-step analyses for the entire list of chalimus-responsive DEPs. For these analyses, we used a list of 1014 DEPs constituted from the union of probes with FC ≥ |2|from ADJ versus PRE and ATT versus PRE lists, and all DEPs from ATT vs. ADJ list ( Figure 1B; Supplementary File S3B). Among 21 level-2 GO BP terms, 11 were found to be enriched (e.g., response to stimulus (GO:0050896), signaling (GO:0023052), immune system process (GO:0002376), and biological adhesion (GO:0022610); Figure 3A). Five GO terms were enriched among 13 level-3 GO terms associated with MF ( Figure 3B), in which many of them were related to 'binding' to different molecules (e.g., ion (GO:0043167) and protein (GO:0005515)), and oxidoreductase activity (GO:0016491). There were 3 level-2 GO CC terms that were enriched out of a total of 13 terms (e.g., 'extracellular region' (GO:0005576; Figure 3C)). It is worth noting that the majority of the GO (>90%) terms were over-represented (Supplementary File S3).  Detailed profiles of GO term distribution of BP and CC domains and their enrichment data for the different segments (i.e., 1, 2, and 3) of DEPs ( Figure 2; FC ≥ |2|) by Blast2GO analysis are available in Supplementary Files S4-S6. Analysis of segment 4 yielded no enriched GO terms. Based on these profiles, multiple GO BP terms associated with stress and/or stimulus-response, and several immune processes were over-represented. Under the CC domain, the majority of the enriched terms were associated with extracellular matrix/region. Several GO MF terms related to O 2 transport (e.g., heme/O 2 binding) were also present in the enriched GO list (Supplementary File S4).

GO/Pathway Term Networks Associated with DEGs
The ClueGO tool [32] was used to perform an integrated GO/pathway term network analysis by functionally organizing the enriched GO (BP, MF, and CC) and Reactome pathway terms. All four segments of DEGs ( Figure 2) were separately analyzed for their up-and down-regulated transcripts, and the results are illustrated in Figures 4 and 5, and Figures S2-S4. Detailed results of these analyses, including information about statistical parameters, clusters, and associated genes, are provided in Supplementary Files S7-S14.
Pathway enrichment results for up-and down-regulated DEGs in segment 1 ( Figure 2; DEGs shared by ADJ versus PRE and ATT versus PRE) are provided in Figure 4 and Supplementary Files S7 and S8. DEGs representing upregulated transcripts of segment 1 formed four main clusters associated with 'cytoskeleton', 'cell cycle and chromosome', 'nucleus', and 'metabolic process' ( Figure 4A). In contrast, GO pathway terms of downregulated transcripts in segment 1 were composed of seven clusters associated with 'extracellular matrix (ECM)', 'cell migration and motility', 'development and differentiation', 'immune system, cell-mediated immunity and defense by lytic activity', 'endoplasmic reticulum (ER), extracellular-related and vesicles', 'phagocytosis', and 'metabolism' ( Figure 4B).   (B) Downregulated DEGs (1277 DEPs). GO enrichment analysis was performed by the ClueGO [32] plugin in Cytoscape [33]. Two databases were used, including GO (BP, MF, and CC) and Reactome pathways for retrieving associated terms. Only the networks and pathways with p < 0.05 are illustrated. Functionally grouped networks with terms as nodes linked by edges based on their kappa score level (≥0.4) are shown. Related-GO terms are grouped and illustrated with distinct color and labeled with the same color. The node size represents the significance of term enrichment. Functionally related groups partially overlap. When a particular GO term is shared by two or more different GO cluster groups, the node is shown by multiple colors. The shape of the nodes indicates the source of the database from where a term was retrieved (ellipse, GO_BP; rectangle, GO_MF; triangle, GO_CC; Reactome, hexagon). The thickness of edges indicates the kappa score (strength of intra-connectivity between cluster groups). Related clusters are shown together within a dotted border and labeled with broad themes for discussion purposes. Refer to the Supplementary Files S7 and S8 for details and high-resolution images. Segment 2 ( Figure 2; DEGs in ATT versus ADJ list) presented six enriched GO term clusters associated with 'cell-matrix and leukocyte migration', 'inflammatory and defense response, acute phase response (APR), and cytokine signaling', 'collagen and ECM degradation', 'wound healing', 'cell proliferation', and 'keratin', among which, the first four clusters were interconnected ( Figure 5 and Supplementary File S9). 'Activation of matrix metalloproteinases' (R-HSA:1592389; corrected p-value 1.86 × 10 −4 ) was the most significant term and possessed the highest percentage of associated genes (19.05%) in this segment.
Enriched GO/pathway terms for segment 3 ( Figure 2) that represented unique DEGs of ADJ versus PRE are shown in Supplementary Figure S2 and Supplementary Files S10, S11. Four major clusters were formed by transcripts with increased expression (Supplementary Figure S2A) and GO terms in these clusters were associated with 'metabolism and gene expression', 'intracellular/nuclear localization', 'chromosome/chromatin', and 'macromolecular metabolism'. The downregulated transcripts of segment 3 were related to 'signaling in immune system', 'cell-mediated and innate immunity', 'vesicle, vacuole and ER', and 'ECM' (Supplementary Figure S2B). Segment 4 ( Figure 2) representing exclusive DEGs of ATT versus PRE mainly featured a large cluster possessing terms of fundamental cellular processes associated with 'chromosome and cell cycle' for upregulated transcripts (Supplementary Figure S3A; Supplementary File S12), and only a platelet-related GO cluster for downregulated transcripts (Supplementary Figure S3B; Supplementary File S13).
Enrichment analyses performed on DEGs with FC ≥ |2| from the entire chalimus-responsive DEPs (n = 1014) with ClueGO are available in Supplementary Figure S4A and Supplementary File S14. 'Regulation of immune system process' (GO:0002682) was the hub for 3 clusters: 'cell-mediated immunity and defense', 'cell migration and adhesion', and 'development and morphogenesis'. Other enriched clusters in this pathway network included 'ECM/collagen degradation', 'metabolism', 'extracellular', 'ER', 'complement system', and 'signaling'. A total of 163 DEPs contributing to the hub GO term, 'regulation of immune system process' (GO:0002682), were used in constructing a heatmap (Supplementary Figure S4B) that revealed the transcript profile of many immune-relevant transcripts (110 DEGs). An overview of GO/pathway term networks analyses demonstrated an overall pattern of immune suppression at both ADJ and ATT sites (by means of downregulating major defense-relevant GO pathways compared to PRE; Figure 4B and Supplementary Figure S2B) and induction of anti-lice responses at ATT sites (by means of upregulating some immune-related GO pathways compared to ADJ; Figures 5 and 6).

Identification of Dysregulated Pathways
The GO/pathway term network analysis identified 249 enriched GO terms. We were able to annotate over half (n = 134) of these enriched GO/pathway terms based on the C5 collection of Molecular Signatures Database (MSigDB; e.g., 'defense response', 'metalloendopeptidase activity', and 'innate immune response'). These 134 enriched pathways and the expression data of contributing

Identification of Dysregulated Pathways
The GO/pathway term network analysis identified 249 enriched GO terms. We were able to annotate over half (n = 134) of these enriched GO/pathway terms based on the C5 collection of Molecular Signatures Database (MSigDB; e.g., 'defense response', 'metalloendopeptidase activity', and 'innate immune response'). These 134 enriched pathways and the expression data of contributing DEGs were analyzed using the Pathifier package [34] to calculate the Pathway Deregulation Score (PDS) for each sample, and samples were clustered with pathways (Supplementary Figure S5A). With respect to PRE control, the degree of dysregulation was demonstrated to be different among ATT and ADJ fin sites; however, the latter two groups were found to have no distinct separation. It was evident that the ADJ fin sites show a higher degree of dysregulation when compared with ATT fin sites (p < 0.05), at least in part, based on the available PDS of 134 GO/pathways, including many immune-relevant GO terms, such as, 'innate immune response' (GO:0045087), 'response to wounding' (GO:0009611), and 'chemokine-mediated signaling pathway' (GO:0070098) (Supplementary Figure S5B-J).

Chalimus-Responsive, Differentially Expressed Transcripts among Different Fin Sample Groups
Heatmaps and clustering analyses provided further visual insights into individual DEPs/DEGs, and their overall and specific expression patterns across sample types. Different DEP sets were subjected to cluster analyses and heatmaps were generated to evaluate if the sample types could be distinguished based on the transcriptional expression profiles ( Figure 6 and Supplementary Figure S6). Direct comparison of ATT and ADJ (90 DEPs) showed that many immune-relevant transcripts were significantly upregulated in ATT sites (e.g., several matrix metalloproteinases (mmps), and their inhibitory counterpart (timp2), cathelicidin (camp), leukocyte cell-derived chemotaxin-2 (lect2), C-X-C chemokine 2 (cxcl2-b), and interleukin 1-beta (il1b); Figure 6) compared with ADJ.
Selected transcripts with important biological roles (based on the enriched functional classes) that demonstrated distinct transcription in terms of pattern and/or magnitude are shown under different functional categories in Tables 1 and 2. Transcriptionally modulated genes were mainly associated with gene expression (transcription and translation), metabolic processes (biosynthesis and catabolism of macromolecules), melanin biosynthesis, ECM organization and disassembly, oxygen transport, redox homeostasis (Table 1), and several components of the immune system ( Table 2). The majority of the selected DEGs were present in segments 1 and 2 ( Figure 2).
Transcription of a large number of homeobox protein family members (hox) was affected, where some of them were upregulated (e.g., hoxa5, hoxc9) while others were downregulated (e.g., hoxa9, hoxc11, hoxd12) in both ADJ and ATT sites compared to PRE (Table 1). CCAAT/enhancer-binding protein beta (cebpb) was induced in ATT sites compared to ADJ sites. Two other transcriptional regulators (i.e., zinc finger protein gfi-1b, gfi1, and peroxisome proliferator-activated receptor delta, ppard) were among the most-induced transcripts in chalimus-attachment sites compared to PRE (Table 1). Sugar or energy metabolism was one of the many metabolic pathways affected during infection, and key enzymes associated with energy metabolism appeared to be upregulated in ATT and ADJ compared to PRE including transketolase (tkt), glucose-6-phosphatase (g6pc), fructose-bisphosphate aldolase C (aldoc), NADH dehydrogenase subunit 1 (nd1), and methyltransferase-like protein 17 (mettl17) ( Table 1). We also found that members of the lipoxygenases (LOX; e.g., aloxe3 (Table 1); alox12, alox15b (Supplementary File S1)) family were upregulated in ADJ and ATT sites compared to PRE control. Three key enzymes involved in melanin biosynthesis (i.e., tyrp1, dct, tyr) were markedly downregulated in both ADJ and ATT sites compared to PRE (Table 1). Substantial changes in the expression of genes encoding structural components of the cytoskeleton and ECM constituents were also observed ( Table 1; Supplementary File S1). Genes encoding collagens (e.g., col10a1, col6a1, col15a1), actins (e.g., tpm1), tubulins (e.g., tuba4a, tubb), fibronectin (fn1), biglycan (bgn), and decorin (dcn) were mostly downregulated in both ATT and ADJ compared to the PRE control, whereas some laminin subunit genes were upregulated in ADJ compared to PRE (lamb1, lamb4) ( Table 1; Supplementary File S1). We observed decreased transcription of regulators of matrix remodeling (tfpi2; Supplementary File S1), and enzymes degrading structural components of ECM (e.g., mmps: 2, 9, 13, 19; serine proteases: htra1, htra3 and cathepsins: cts) in ATT and ADJ compared to PRE (segment 1; Table 1). However, it should be noted that both mmps (mmp9-a, mmp13-b, mmp13-c, mmp14) and timp2-b demonstrated an induced expression in ATT compared to ADJ (segment 2). Some genes associated with oxygen transport (hemoglobin subunits: hba, hba4, hbb) and iron metabolism (5-aminolevulinate synthase, alas2) also presented reduced transcript levels in lice-infected salmon (ATT and ADJ) compared to PRE. However, an iron-binding protein of this group, cytoglobin-2 (cygb), was upregulated in ADJ sites compared with PRE (segment 3; Table 1). Several members of the antioxidant enzyme family were found to be affected during lice infection. The transcription of thioredoxin-related (txn-b, txndc5), peroxiredoxin (prxl2a), glutathione peroxidase (gpx7), glutathione S-transferase (gsta), and ferritin (fth1, frim) genes was decreased in ADJ sites compared to PRE. Although the transcription of many of these genes in ATT compared with PRE was also reduced (segment 1; txndc5, prxl2a, gpx7, gsta), expression of txn-b was higher in the ATT site compared to the ADJ site (segment 2; Table 1).
In the present work, dramatic changes in various immune-relevant pathways and mechanisms were seen post-lice exposure ( Table 2). Decreased mRNA abundance of various pattern recognition receptors (PRRs) was observed (different cluster differentiation antigens, cd93, cd209, cd302; galectin-1, lgals1; proteins involved in mannose recognition, mbl2, and mrc1; two C-type lectin members, clec4m and clec4e). Chalimus-attachment induced the expression of lect2, C-X-C chemokines (cxcl2, cxcl11), receptors (e.g., cxcr1), and interleukins (il1b and il11) in ATT compared to ADJ. A C-C chemokine (ccl4) and a receptor (cmklr1) were downregulated during infection (both in ATT and ADJ) compared to PRE. Genes encoding different complement factors (cf such as cfh, c4a, cfd, but not cfb) also showed reduced expression in ATT and ADJ compared to PRE (segment 1; Table 2). Further, downregulation was noticed for complement receptor, c3ar1, along with cd59, in ATT and ADJ compared to PRE. Players in the coagulation cascade demonstrated a differential response to sea lice. Coagulation factor X (f10) was induced; however, f3 was downregulated in ATT and ADJ compared to PRE (segment 1; Table 1). Meanwhile, increased transcription of f5 and plasminogen activator inhibitor 1 (serpine1) at ATT sites was observed compared to ADJ (Table 1). A panel of mRNAs encoding proteins involved in inflammation and APR (e.g., vascular cell adhesion protein 1, vcam1; midkine, mdk; C-reactive protein, crp; prostaglandin E synthase, ptges) was found to be downregulated in ATT and ADJ compared to PRE. However, the gene encoding prostaglandin D2 receptor 2 (gpr44) was upregulated in ADJ versus PRE fin samples ( Table 2). The downregulation of viral-responsive genes (e.g., rsad2-a (alias viperin); and two interferon-induced proteins, ifi44 and ifit5, transcripts encoding TRIM family members, trim8 and trim58) in both ATT and ADJ compared to PRE, was a characteristic response to chalimus infection. Nevertheless, interferon α3 (ifna-a) mRNA increased in ATT and ADJ compared to PRE ( Table 2). A set of innate immune genes demonstrated significantly reduced transcription either in the ADJ site alone (complement component 1 Q binding protein, c1qbp and camp-b) or both in ADJ and ATT sites (putative defense protein HDD11, hdd11; CD99 antigen-like protein 2, cd99; complement c1q tumor necrosis factor-related protein 6, c1qtnf6; fibrinogen-like protein 1, fgl1; lipopolysaccharide-binding protein, lbp; lipopolysaccharide-induced tumor necrosis factor-alpha, litaf ) compared to PRE ( Table 2). In contrast, NLR family CARD domain-containing protein 3 (nlrc3) and complement C1q tumor necrosis factor-related protein 3 (c1qtnf3) showed induced expression in ATT and ADJ compared to PRE. Two camp paralogs, fibroblast growth factor 1 (fgf1) and arginase-2 (arg2) genes, demonstrated significantly higher expression in ATT to that of ADJ. Genes encoding important components of adaptive immunity, such as two major histocompatibility (MH) class I antigens (hlab and hlah) and an immunoglobulin chain (e.g., iglc3), were also significantly downregulated in ATT and ADJ compared to PRE. itga10 Integrin alpha-10 −2.6 NS NS 3 ctgf Connective tissue growth factor (CCN family member 2) Hemoglobin subunit alpha (and -4) Coagulation factor III (tissue factor) −3.9 −2.3 NS 1 (Some of these also take part in acute phase response (APR)) Tissue factor pathway inhibitor 2 precursor

QPCR Confirmation of Selected Transcripts
A subset of chalimus-responsive transcripts (37 microarray-identified and one additional MMP family member) was chosen for confirming the microarray results by QPCR. These transcripts were associated with five different functional themes including (1) ECM degradation, tissue repair/remodeling, and wound healing, (2) immunity and defense, (3) chemotaxis and signaling, (4) antiviral response, and (5) antioxidant activity and redox homeostasis (Table 3). In addition to the samples used in microarray, we included 4 more samples/group in the QPCR analysis (n = 10 in total/group). Log 2 -transformed microarray and QPCR FC showed a significant positive correlation (R 2 = 0.76, p < 0.0001; Figure 7A). Individual relative quantity (RQ) profiles for QPCR with FC values from both QPCR and microarray for each theme are presented in Figures 7-10. Based on inset tables, the FCs obtained from two techniques presented some divergence across genes of interest (GOIs), but very similar patterns in the direction of change (i.e., up-or down-regulation), in general (Figures 7-10). Specifically, 32 transcripts (86.5%, except bgn-a, c1qtnf3, cebpb, ifna-a, and txn-b) were confirmed as significantly (p < 0.05) differentially expressed in at least one comparison among three sample groups (i.e., ADJ versus PRE, ATT versus PRE, or ATT versus ADJ), from which they were identified as informative DEPs by SAM (Supplementary File S1).
QPCR profiles of transcripts involved in ECM degradation, wound healing, and tissue repair/remodeling (theme 1) are shown in Figure 7. The majority of these transcripts demonstrated lower expression in infected animals compared to PRE control animals. Significant differences (p < 0.05) were observed in both ATT and ADJ fins for timp2-b, timp2-c, mmp13-d, mmp2, mmp19, mmp9-a, mmp9-b, and bgn-b. Except for mmp9-a, mmp20, and bgn paralogs, the expression of other transcripts in ATT was higher compared to ADJ, although this was only significant (p < 0.05) for mmp13-c based on QPCR results ( Figure 7E). It was evident that, in general, transcripts in theme 1 were expressed at a lower level in ADJ compared to PRE (and ATT for some transcripts) (Supplementary Figure S7A).
Conversely, the transcripts of theme 2 that are immune/defense (not including antiviral response)-relevant demonstrated two distinct patterns of expression (Figure 8), as revealed by a heatmap (Supplementary Figure S7B). A subset of transcripts (i.e., camp-a, f5, and c1ql2) were more highly expressed in ATT compared to ADJ and PRE. Significant differences (p < 0.05) in mRNA levels between ATT and ADJ were only detected for camp-a and f5 by both QPCR and microarray ( Figure 8A,B). Although it showed a similar trend with the microarray data, clqtnf3 transcription did not differ significantly among sample groups ( Figure 8C). The mRNA level of c1ql2 was significantly higher in ATT compared with PRE (p < 0.05; FC > 9; Figure 8D). The second set, comprising camp-b, lgals1, cd59, cd83, crp, and hdd11 transcripts, showed significant (p < 0.05) lower expression in ADJ and/or ATT compared with PRE ( Figure 8E-J). QPCR failed to detect any significant difference in the expression of cebpb ( Figure 8K).
Many transcripts of both themes 3 (chemotaxis and signaling) and 4 (antiviral response) showed similarly repressed expression in ADJ and ATT, as of theme 1, when compared with PRE samples, except for cxcl2-b, ccrl1 and ifna-a ( Figure 9, Supplementary Figure S7C). Chalimus-attachment significantly (p < 0.05) upregulated the transcript level of cxcl2-b in ATT compared to ADJ ( Figure 9A). While the expression level of lect2-a, ccl17, cxcl11 and rsad2-a was significantly (p < 0.05) downregulated in ADJ and ATT, compared with PRE ( Figure 9B-D,F), ccrl1 was significantly (p < 0.05) upregulated ( Figure 9E) in both ATT and ADJ compared to PRE controls. According to QPCR data, ifna-a was transcriptionally non-responsive to chalimus infection ( Figure 9G). Table 3. Summary of details about the real-time quantitative polymerase chain reaction (QPCR)-analyzed genes of interest (GOIs) from different functional themes.  Figure S2B fibronectin fn1 ECM assembly, cell motility, adhesion, wound healing Extracellular space, ECM organization, Acute-phase response 1, 2 Figure 4B, Figure 5 biglycan-like bgn Collagen fiber assembly, ECM organization, matrix mineralization Extracellular space, ECM organization, Collagen-containing ECM 1 Figure 4B 2. Immunity and defense (Not including antiviral response) cathelicidin camp Antibacterial activity Defense response to bacterium, Antimicrobial humoral immune response mediated by antimicrobial peptide 2 Figure 5 beta-galactoside-binding lectin lgals1 Regulating apoptosis, cell proliferation and cell differentiation Extracellular space, Regulation of immune system process 1 Figure 4B CD59 glycoprotein cd59 Complement MAC inhibition Regulation of immune system process, Cell surface receptor signaling pathway 1 Figure 4B CD83 cd83 Antigen presentation and lymphocyte activation Regulation of immune system process, Leukocyte activation, Positive regulation of immune system process 1 Figure 4B complement C1q tumor necrosis factor-related protein 3 c1qtnf3 Regulatory role in immune system $ Negative regulation of inflammatory response, Positive regulation of cytokine secretion, Immune system process 1 Figure 4A complement C1q-like protein 2 c1ql2 Regulate excitatory synapses $ Extracellular region, Immune system process 4 Supp. Figure S3A C-reactive protein crp Inflammation, acute phase response, defense Leukocyte mediated immunity, Regulation of immune system process, Defense response 1 Figure   Chemotactic activity, chondrocyte proliferation $ Response to stimulus, Immune system process 2 Figure 5 C-C motif chemokine 17-like ccl17

T-lymphocytes trafficking and activation
Response to cytokine, Cellular response to cytokine stimulus 3 Figure 2B chemokine receptor 4-like ccrl1 Regulator of chemotaxis $ Response to stimulus, Immune system process 1 Figure 4B C-X-C motif chemokine 11-like cxcl11 Chemotactic for activated T-cells, neutrophils or monocytes Myeloid leukocyte migration, Neutrophil chemotaxis, Chemokine-mediated signaling pathway 2 Figure 5 4. Antiviral response radical S-adenosyl methionine domain-containing protein 2, alias (viperin) rsad2 Antiviral Immune system, Regulation of immune system process 1 Figure 4B interferon a3 ifna * Antiviral Immune system process, Response to stimulus 1 Figure  NADPH oxidase, redox homeostasis $ Superoxide-generating NADPH oxidase activator activity, Respiratory burst 1 Figure 4B 1 Based on functional classes, candidate genes are categorized under five themes for discussion purposes. 2 Name(s) or alias(es) obtained from annotation. Refer the Supplementary File S1 for additional details. 3 Official gene symbols (version 1) are based on multiple annotations and the majority of them are represented in HUGO Gene Nomenclature Committee (HGNC) (https://www.genenames.org/) and/or GeneCards (https://www.genecards.org/) databases. 4 General overview of the function of each candidate gene extracted from the UniProtKB database (https://www.uniprot.org/uniprot/). 5 Selected enriched GO terms associated with a GOI obtained from enrichment analyses. Additional details are available in Supplementary files. 6 In which segment (refer to Figure 2) of the microarray data list, the majority of the DEPs representing this transcript is present. # This was not identified by the microarray in this study. * Putative gene symbol of ifna (LOC106607463). $ One or more terms were obtained from UniProtKB/Blast2GO analyses.  Finally, QPCR (and microarray) results showed that chalimus infection robustly and significantly (p < 0.05, except txn-b, which showed non-significant QPCR results) suppressed the expression of all the studied transcripts representing antioxidant and redox system in theme 5 ( Figure 10). Both QPCR and microarray did not find significant (p < 0.05) differences between ATT and ADJ sites in fins for gsta, gpx7, prxl2a, prdx1-a, and ncf2. However, as illustrated in the heatmap (Supplementary Figure  S7D), transcription of theme 5 genes appeared to be more highly suppressed in ADJ compared to that of ATT (with respect to PRE; see FCs).

Multivariate and Correlation Analyses Based on Gene Expression and Phenotype Data
To examine the similarities and dissimilarities among the fin sample groups, and to understand the relationship among different phenotypic parameters, we performed a PCoA. As shown in Figure 11A, PCoA explained 63.6% of the variation among sample types. While PCO1, driven by many transcripts (e.g., mmp9-a, mmp13-d, cd59, ncf2, rsad2, and prdx1a) segregated PRE control salmon fins from infected salmon fins, PCO2 was greatly influenced by f5 and camp-a, and separated fin samples of infected animals based on sampling sites (i.e., ATT versus ADJ). Permutational multivariate analysis of variance (PERMANOVA) indicated that the sample groups are significantly different from each other (p = 0.0001). Higher expression of the majority of the QPCR-assayed genes in the PRE group was also evident from PCoA. To complement PCoA, we conducted Pearson correlation analyses. Seven different phenotypic parameters and gene expression data were used in this correlation analysis, and the results are tabulated in Figure 11B. Chalimus load negatively correlated with many genes that were QPCR-tested (e.g., ifna-a, cd59, crp, ccl17, txn-b, prxl2a, and prdx1a). Conversely, cd83 revealed a positive correlation with the chalimus load. Interestingly, viscerosomatic index (VSI) was positively and significantly (p < 0.05) correlated with 22 genes tested, while negatively and significantly (p < 0.05) correlated with ccrl1 ( Figure 11B). Figure 11. Principal coordinate analysis (PCoA) and correlation analysis of fin transcript expression data from QPCR with phenotypic factors. Standardized and log 2 -transformed data, respectively, were used in PCoA and correlation analyses. (A) Similarities and dissimilarities among treatment groups revealed by PCoA performed by PRIMER 6.1.15. (B) log 2 -transformed data were subjected to factor analysis using IBM SPSS Statistics (v25). Matrix represents the Pearson correlation coefficients between 8 phenotypic factors (columns; e.g., chalimus count) and both fin gene expression and the phenotypic factors (rows). Positive and negative correlations are indicated with green and red, respectively. A number for the corresponding correlation is displayed if the correlation is significant (p < 0.05). TCC, total chalimus count; W, weight; L, length; HSI, hepatosomatic index (HSI = 100×(liver weight/W)); SSI, spleen-somatic index (SSI = 1000× (spleen weight/W)); VSI, viscerosomatic index (VSI = 100× (viscera weight/W)), and K, Fulton's condition factor (K = 100× (W/L 3 )).

Discussion
Anti-parasitic responses in teleosts have received less attention compared to studies focused on responses to bacterial and viral infections [35]. There is little published information available on early immune responses of Atlantic salmon against sea lice (L. salmonis) [24,25]. We used a moderate-density infection model (ca. 0.15 L. salmonis/g fish) in our study. Johnson and Albright [36] reported an abundant presence of early life stages (<5 dpi) of L. salmonis on gills of Atlantic salmon, while later life stages (10-20 dpi) preferentially settled on fins when compared to gills and other body surfaces. A comparative study reported that L. salmonis and C. elongatus chalimi were abundantly attached to fins compared to body and head, with the lowest abundance in gills under natural conditions in Atlantic salmon, whereas gills were a major attachment location in experimental models compared to fins and body [37]. Chalimi distribution in our infection trial was in agreement with these studies [36,37], and a considerable number of chalimi settling on the gills could be an artifact of laboratory conditions [37]. By being the initial and preferential settlement sites of chalimi, fins were targeted for gene expression profiling in the present study.
Some studies have examined the spectrum of relative susceptibility in different salmonid species to L. salmonis infection. Results of these studies highlighted that some salmonids (e.g., chinook (O. tshawytscha), chum (O. keta) or sockeye (O. nerka) salmon) in addition to the Atlantic salmon are highly susceptible, whereas others are more resistant (e.g., coho and pink salmon) [8,9,36,[38][39][40]. The high infection rate and prolonged retention of parasites throughout the trial period in our study (i.e., chalimi versus pre-adults counts) confirmed Atlantic salmon's inability to reject sea lice and its status as a susceptible host [2].
It is difficult to determine the gross changes in whole tissue transcriptome in response to parasites using the conventional targeted QPCR approach. However, microarray-based high-throughput platforms are used to screen tens of thousands of genes and hundreds of pathways simultaneously and to identify and characterize known and putative markers. Herein, we focused on 8 dpi ( Figure 12A), at which time copepodids had developed to the chalimus I stage. We employed microarray techniques to identify the differentially expressed transcripts and elucidate the pathways and mechanisms involved in early response against sea lice attachment on Atlantic salmon fins.
The results of two independent approaches (microarray and QPCR) were generally in significant concordance. However, we observed some transcripts that were differentially expressed as per SAM (e.g., bgn-a, c1qtnf3, cebpb, ifna-a, and txn-b) but with no significant differences in RQs from QPCR across treatment groups, or vice versa (Figures 7-10). This divergence may have been caused by different statistical approaches used in microarray (i.e., SAM (a modified t-test) for pairwise comparison of treatment groups) and QPCR (i.e., analysis of variance (ANOVA) compared the means of three independent groups). The increased number of biological replicates enrolled in the QPCR study could also be a source of this disagreement.
Transcriptional differences between paralogs of bgn, camp, mmp9, mmp13, and timp2 were addressed in the current study. Paralogs of some genes demonstrated similar mRNA expression profiles (e.g., timp2, mmp13 (paralogs b and d), and bgn), whereas others presented distinct transcript expression profiles (e.g., mmp9, camp) (Figures 7-10). The entire teleost lineage has undergone three rounds (1R, 2R, 3R) of whole-genome duplication (WGD) that shaped their evolution. While 3R was a teleost-specific WGD, a salmonid-specific fourth WGD (4R) occurred in the common ancestor from which different sub-lineages of salmonids evolved [41,42]. These WGDs have generated a complex and diverse repertoire of paralogs from ancestral genes in salmonids. The potential fates of duplicated genes post-WGD include subfunctionalization (each paralog retaining a subset of its original ancestral function) or neofunctionalization (a paralog acquiring a new function) [41,42]. The paralog-specific transcript expression patterns from our study provide evidence for gene duplication and divergence. However, whether these paralogs were subfunctionlized or neofunctionalized is not known and requires functional studies with the encoded proteins.
Previous research found that local response in fin/skin to lice was more diverse and complex when compared with the systemic response (such as in spleen) [24]. In the current study, 6568 DEPs were identified ( Figure 1B; representing 3571 DEGs), and the larger fraction of DEPs was present in ADJ versus PRE compared to ATT versus PRE and ATT versus ADJ. Quantitative estimation of pathway dysregulation indicated that the ADJ group presents higher PDS compared to ATT and PRE groups for the majority of the GO terms, including many immune-relevant GO terms (Supplementary Figure S5), and this is in agreement with the relative distribution of DEPs among sample groups ( Figure 1B), PCoA of DEPs ( Figure 1C), and the RQ values of the majority of GOIs determined by QPCR (Figures 7-10). Several aspects of gene expression profiles observed in our study provided evidence suggesting that specific immune signaling pathways were modulated in infected fins (to reject sea lice). It was interesting to note that~30% of the entire set of chalimus-responsive DEPs (1928, segment 1; Figure 2) were common among ATT and ADJ compared to PRE, and each of them varied only in the level of mRNA expression, not in their direction of modulation (Supplementary File S1).
Our expression profiling and enriched GO/pathway term analyses collectively uncovered that chalimus infection caused transcriptional changes in an array of biological processes (Tables 1 and 2), some of which are discussed below.

Chalimus Infection Influences Transcription Machinery, Energy Metabolism, and Melanin Biosynthesis
As noticed in our GO enrichment analyses (gene expression (GO:0010467); Supplementary File S7), a wide range of homeobox transcription factors (TFs) has been seen to be modulated by stress conditions, including infection by sea lice in salmon [25,26]. Although members of the HOX family were originally found to determine cell fate and regulate organ development [43], recent discoveries have implicated them in wound healing and tissue remodeling [44]. CEBP is a family of TFs governing immune and inflammatory responses [45], and cebpb (encoding CEBPβ protein) was proposed as a lice-responsive marker in salmonid skin by many studies [23,39]. This transcript was upregulated in skin of Atlantic salmon as early as 24 hpi with adult L. salmonis [38]. NF-κB is another vital synergetic TF of CEBPβ that activates inflammatory responses, and nfkb1 and nfkb2 mRNAs were significantly upregulated in ADJ sites compared to PRE (~1.5-fold; Supplementary File S1). The upregulation of gfi1b in ATT and ADJ compared to PRE may lead to transcriptional repression in erythroid cells, as shown in mouse [46]. PPARD is a member of the ligand-inducible PPAR TF family whose increased transcription might enhance fatty acid oxidation and energy uncoupling in muscle and adipose tissues, and suppress macrophage-derived inflammation [47,48]. Modulated expression of transcripts encoding all these TFs, which are involved in a multitude of functions, such as cytokine production, inflammation, wound healing, and hematopoiesis, was suggestive of major changes in multiple downstream events in our study.
Based on our GO/pathway term analyses, enriched GO term clusters represented by the upregulated DEGs in the ATT and ADJ compared to PRE (e.g., 'metabolism', 'cell cycle and chromosomes', and 'cytoskeleton') implied that fin regeneration processes might be in progress in lice-infected fins ( Figure 4A, Supplementary File S7). Gross changes in metabolic processes of various macromolecules (sugars, amino acids, and nucleic acids) were apparent, particularly in ADJ and ATT sites compared to PRE (Group 3; e.g., heterocycle metabolic process (GO:0046483), Supplementary File S7). Induced transcription of genes involved in sugar and lipid metabolism in infected animals compared to PRE suggested an increased energy mobilization during infection. Energy metabolism-related genes were also present among DEGs in Caligus-infected salmon fins [25]. Upregulated expression of several transcripts encoding LOX enzymes (particularly, aloxe3 and alox12; Supplementary File S1) that play key roles in biosynthesis of oxylipin signaling molecules [49] suggested that the salmon attempted to modulate the permeability barrier function of fin against sea lice.
Our enrichment analyses indicated that melanin-related pathways are downregulated in lice-infected salmon (e.g., melanosome (GO:0042470); Figure 4B, Supplementary File S8). Kittilsen et al. [50] showed that hyper-melanized salmonids demonstrate a reduced physiological and behavioral response to stress. In addition, the same group later reported that melanin-based pigmentation significantly correlated with the immunocompetence of Atlantic salmon and its resistance to sea lice [51]. A recent study found hyperpigmentation and migrating pigment bodies at the wound surface and implicated it with wound repairing processes in Atlantic salmon [52]. It may be argued that parasite-mediated inhibition of melanin synthesis could be a potential mechanism supporting the lice pathology.

Sea Lice Modulate ECM Integrity and Catabolism, Wound Repair, Inflammation, Acute Phase Response, and Coagulation
Complex dynamics in transcriptional expression of different structural and functional components of ECM were evident from our GO/pathway term and expression analyses (e.g., group 5 representing ECM; Supplementary File S8). Notably, transcripts encoding several structural elements of ECM, such as fibronectin, collagen chains, laminin, and proteoglycans (BGN, DCN), were downregulated, suggesting impaired tissue repair and wound healing cascades at the infection site. Beyond its structural roles in ECM, BGN orchestrates signaling networks involved in inflammation and immunity [53]. Meanwhile, DCN is a multifunctional protein involved in anti-fibrosis by interacting mainly with TGFβ and other matrix molecules and cytokines (Reference [54] and references therein). Both BGN and DCN promote cell migration and are involved in the remodeling process [55]. Downregulation of fn1 (fibronectin) in liver throughout the lice-infection led Skugor et al. [23] to suggest an impaired ability to heal wounds in salmon. Moreover, mRNAs encoding several lysosomal (serine (HTRA1, HTRA3) and cysteine (cathepsins)) proteases that are involved in ECM catabolism were also downregulated in ATT and ADJ compared to PRE. Annexins are important ECM proteins that play multiple roles in membrane scaffolding, cell division, signaling, and apoptosis [56], and genes of several annexin family members (e.g., anxa13) were downregulated in ATT and ADJ compared to PRE. As a collagen-specific chaperone, SERPINH1 plays roles in collagen biosynthesis and structural organization [57]. serpinh1 mRNA demonstrated decreased abundance in ATT and ADJ compared to PRE.
The TGFβ proteins (a family of growth factors and cytokines), their receptor complexes, and associated proteins together build an intracellular transmission cascade that governs proliferation, differentiation, and migration of cells and determines the healing dynamics [58]. As noticed by Krasnov and colleagues [26], chalimus-attachment appeared to inhibit the TGFβ signaling network based on suppressed transcription of ligand (tgfb3), a receptor (tgfbr2), and several regulatory proteins (lrrc32, tgfbi [15], and ltbp3) of this signaling pathway (Supplementary File S1). Transcription of other growth factors and regulators that operate in conjunction with TGFβ, namely CTGF, FGF1 [59], and SERPINF1 [26], was also impacted in our study. Overall, these transcriptional responses may result in a lethargic wound repairing process.
There is an equilibrium in expression and activity between MMPs and their endogenous inhibitors (TIMPs) in the ECM microenvironment [60]. This balance primarily determines the dynamics of ECM degradation and remodeling. We identified a set of differentially transcribed mmps (collagenases and gelatinases) along with two of their inhibitors (timp2 and timp3). Paralog-specific QPCR assays were developed, and results indicated distinct expression patterns (Table 1, Figure 7). Although only QPCR confirmed for mmp13-c, the microarray results suggested a more pronounced expression of several mmps (e.g., mmp9-a and mmp13-b) and timp2-b in ATT compared to ADJ. Previous evidence from different groups suggests that lice-infection chronically damages the tissues, resulting in a hallmark aberrant expression of mmps [8,23,38,39]. Two enriched GO terms related to wound healing (i.e., response to wounding (GO:0009611) and wound healing (GO:0042060)) represented many associated genes that were downregulated in infected salmon compared with PRE ( Figure 4B, Supplementary File S8). Collectively, these data suggest that the lice infection influences ECM catabolism, remodeling, and wound repair events in Atlantic salmon fin tissue.
Both microarray and QPCR indicated downregulated transcription of the gene encoding CRP, a well-known acute phase protein (APP) [61], in infected fin (Table 2, Figure 8I). Pro-inflammatory cytokines (e.g., IL1β, CXCL8) induce APR through the activation of NF-κB and CEBPβ [62]. Despite the upregulation of il1b mRNA, decreased expression of crp transcript in infected fin compared to PRE might have resulted in a suppressed APR during lice-infection. Limited local expression of crp in lice-attachment sites has been implicated with susceptibility of S. salar to sea lice based on comparison with resistant salmonids [24,38]. Decreased transcription of mdk implied an impaired tissue repair process as MDK enhances inflammation in mammals by promoting the migration of inflammatory leukocytes [63]. We found that transcription of ptges was also downregulated in ADJ and ATT compared to PRE (Table 2). Reduced transcription of genes encoding prostaglandin synthase enzymes (e.g., PTGES) in Atlantic salmon skin during lice-infection has previously been reported [24,38,39]. All these events may plausibly dysregulate inflammation by restricting and/or delaying it, which in turn may result in chronic wounds and higher susceptibility to parasitization [10,64]. This may lead to the cutaneous lesions and loss of osmotic balance in wounded areas, causing a greater amount of distress to infected salmon [65].
Some components of the coagulation cascade also play roles associated with APR (e.g., fibrinogen, plasminogen, and plasminogen activator inhibitor 1 (PAI-1)) [66] and inflammation (e.g., coagulation factor III and thrombin) [67]. As found in previous studies [24,68], members of the coagulation cascade, and hence hemostasis, were found to be affected by lice-infection in our study. We observed the downregulated transcription of f3, the primary initiator of the extrinsic coagulation pathway [69], and tissue factor pathway inhibitor 2 (tfpi2; [70]), in ATT and ADJ compared with PRE ( Table 2). Although the blood-feeding characteristic of sea lice has been documented for matured developmental stages of sea lice [71], it is not clear whether or not the lice larvae possess any mechanisms to maintain a steady flow of host blood at feeding-sites. Conversely, a few genes encoding proteins with regulatory roles in the coagulation pathway, coagulation factor V (co-factor of factor Xa) and PAI-1 (that inhibits fibrinolysis), were induced in ATT compared to ADJ (Table 1). These findings suggest that hemostasis was dysregulated by chalimi, whereas salmon attempted to restore the coagulation mechanism at the chalimus-attachment site, perhaps to prevent the lice from blood-feeding.

Chalimus-Induced Changes in O 2 Transport and Redox Homeostasis
The heme biosynthetic pathway is absent in many parasites [72]. Sea lice are hematophagous species that acquire iron from their hosts for their nutritional requirements. A 'nutritional immunity' concept hypothesized that hosts restrict the access of iron or heme groups by pathogens [73]. Evidence from different studies at the salmon-louse interface provided contrasting results. Heggland et al. [74] recently compared the gene expression of attached chalimi in gills and skin of Atlantic salmon and found elevated transcripts of genes important for absorption, storage, and/or transportation of iron and heme, digestive and detoxification enzymes, and anti-clotting elements in gills. Iron metabolism was shown to be modulated by means of increased transcriptional activities of genes associated with heme biosynthesis (e.g., alas) and iron transport (fth1, frim, haptoglobin (hp), and hepcidin (hamp)) coupled with decreased expression of enzymes degrading the heme (heme oxygenase, hmox1) [75]. However, other authors found coordinated suppression in systemic and/or local mRNA expression of hemoglobins, alas2, and hmox1 [23,26]. In the current study, several probes representing hemoglobin subunits and alas2 demonstrated decreased mRNA abundance post-infection compared to PRE (Table 1). Sustaining a lower blood flow and heme enrichment cannot be ruled out as a potential protective mechanism of the host to deprive the parasite of iron supply [26]. Additionally, proteins encoded by some of these genes (e.g., hp, ferritin, and hamp) are also considered to take part in APR [66]. The elevated transcript level of another oxygen-carrying protein, cytoglobin-2 (cygb), in ADJ compared to PRE (Table 1), suggested the existence of potential alternative mechanisms to protect the host tissues from hypoxia.
Excessive stress conditions, including pathogenic infections, could result in the production of a massive amount of reactive oxygen species (ROS). Although ROS play crucial roles at optimal physiological levels, they could cause damage to cells at higher concentrations [76]. Enzymatic antioxidants neutralize the excessive ROS and restore the cellular redox equilibrium. Many genes encoding antioxidant enzymes were transcriptionally suppressed (e.g., gsta, gpx7, prxl2a, and prdx1-a) by chalimi in infected animals compared to PRE, indicating that redox homeostasis was disturbed by the sea lice infection, as reported earlier [77]. However, in agreement with previous studies in salmonid skin [23,39], our microarray results revealed upregulated txn-b expression in ATT compared with ADJ (Supplementary Figure S4B). Collectively, these results suggest that TXN could be a crucial player in maintaining redox homeostasis during lice-infection.

Impact of Lice Infection on Expression of Immune-Relevant Transcripts
Although an overall downregulation of immune genes was a hallmark of lice pathology ( Figure 4B, Supplementary Figure S4B; e.g., a large number of GO terms in groups 13 and 19 in Supplementary File S8), some immune markers were upregulated in ATT compared to ADJ (Figure 6; e.g., GO terms of groups 7 and 10 in Supplementary File S9).

Pattern Recognition Receptors (PRRs)
PRRs play vital roles in innate immunity by recognizing pathogen-or damage-associated molecular patterns (PAMPs or DAMPs) [78]. These PRRs sense danger signals and initiate complex cascades culminating in the production of effector molecules, such as cytokines, ROS, antimicrobial peptides, growth factors, and complement proteins [79]. In line with previous studies [8,[23][24][25][26], we found an overall transcriptional downregulation of several PRRs in infected salmon compared to PRE (e.g., cd209, cd302, mbl2, and lgals1). Expression of mrc1 that encodes a mannose receptor was downregulated in ADJ compared with PRE ( Table 2). As a member of the calcium-dependent C-type lectin receptor family (CLR) [80], CLEC4E is expressed in macrophages and dendritic cells [81]. The transcript encoding CLEC4E was downregulated in ADJ (compared with PRE) and upregulated in ATT compared with PRE. Two other members of CLRs, CD209 and CD302, are abundantly expressed in dendritic cells, and implicated with phagocytic pathogen-recognition [82] and dendritic cell migration [83], respectively. Mannose-binding lectins (e.g., MBL2 and MRC1) also belong to another subgroup of CTLs [84]. A recent study found a strong correlation between mrc1 transcription and C. rogercresseyi chalimus-load in Atlantic salmon [25]. Galectin-1 (encoded by lgals1) is a glycan-binding protein that acts as a master regulator by controlling proinflammatory cytokine production, neutrophil trafficking, and eosinophil migration. Based on these findings, we speculate that lice-infected salmon could be vulnerable to secondary infections due to the compromised pathogen-sensing ability.

Cytokine/Chemokine Signaling
The present study found altered transcriptional patterns for a large number of cytokines/chemokines and their receptors. A subset of these genes demonstrated downregulation in both ATT and ADJ compared to PRE (e.g., ccl4, cmklr1) (segment 1; e.g., signal transduction (GO:0007165), Supplementary File S8). CCL4 (MIP1β) is a potent lymphocyte chemoattractant, and activation of chemotaxis by CCL4 towards the injury site has been demonstrated in mammals [85,86]. Primarily expressed by dendritic cells and macrophages, CMKLR1 transduces the signals by chemerin, which is an antimicrobial molecule expressed in human skin and induces transmigration of various immune cells [87]. The CCR4-CCL17 receptor-ligand pair showed decreased expression in ADJ compared to PRE (segment 3; Figure 2), indicating coordinated changes in its signaling axis. In mammals, CCR4 is the receptor for CCL17 and is expressed in skin-homing distinct T cell subsets, including activated T cells, Th2 cells, and Treg cells, and implicated in skin-associated immune responses [88].
In contrast, the second group of chemokine-mediated genes showed significantly higher transcript abundance in ATT compared to ADJ (with GO annotation chemokine-mediated signaling pathway (GO:0070098), segment 2 ( Figure 2); Supplementary File S9): CXCR1 coordinates neutrophil trafficking by liganding with CXCL6 and CXCL8 in inflamed tissues [89]. Several LECT2 homologs have been identified from various fish species and demonstrated to have conserved roles as a chemoattractant. LECT2 was identified as an APP during bacterial infection in zebrafish [90] and presented an mRNA upregulation in a lice density-dependent manner in both wild and farmed Atlantic salmon skin at 24-26 dpi [91]. The disagreement in lect2-a mRNA expression in our study and that of Gallardi et al. [91] might be influenced by differences in investigated tissues and sampling time points. Human CXCL2 (MIP2) is chemotactic for polymorphonuclear leukocytes, especially for neutrophils, and transcribed at the sites of wound repair [92]. We noted an interesting expression pattern for cxcl11 (upregulated in segment 2) and its receptor cxcr3 (downregulated in segment 1, Figure 2; Supplementary File S1). In zebrafish, the CXCR3-CXCL11 axis mediates macrophage recruitment against bacterial infection [93]. IL1β is a pivotal proinflammatory cytokine with diverse physiological roles and its major contribution in regulating the inflammatory process is conserved in fish [94]. Despite its transcriptional upregulation in ATT compared to ADJ, the decoy receptor of IL1β (i.e., IL-1RII) was downregulated in ADJ compared to PRE (Supplementary File S1). Mammalian IL11 signaling plays a crucial role in thrombopoiesis [95] and we found local upregulation of il11 in ATT compared to ADJ. Taken together, our data suggest that transmigration and trafficking of different immune cells in salmon fin are dramatically altered through the transcriptional modulation of cytokines and receptors at the site of infection.

Complement Pathway
We documented coordinated suppression of components of the complement cascade (c4a, cfh, and cfd), receptors (cd93, c3ar1, and c5ar1) and some regulators (cd59) of complement activation in this study. As an integral part of the immune system, the complement pathway detects pathogens, alarms the defense system, and eliminates pathogens. Moreover, it serves as a bridge between adaptive and innate immunity and orchestrates acquired immune responses [96]. Receptors of C3a and C5a (potent chemotactic and pro-inflammatory) peptides are expressed by macrophages and shown to participate in signal modulation of the complement system [97]. Mammalian CD93 is a C1q-receptor and implicated with intracellular adhesion, phagocytosis, and inter-cellular interactions [98]. Downregulation of complement-related genes in skin during the lice-infection has been previously reported [24]. Collectively, our findings indicate a potential disruption in the routine function of complement apparatus. Conversely, a member of the C1q family (c1ql2; [99]) was transcriptionally upregulated in infected animals compared to PRE ( Figure 8D), but its association with lice infection is unclear.

Antiviral Responses
In the present microarray study, several viral-induced transcripts (e.g., ifi44, ifit5, trims, and rsad2) were downregulated in ATT and ADJ fin of lice-infected animals compared to PRE (Table 2). In the context of lice-infection, suppressed antiviral responsive pathways were evident in anterior kidney and skin of different salmonid species [8,26,39]. Barker et al. [13] recently demonstrated decreased transcript levels of antiviral-effector proteins (e.g., mx, mhc-Ib, galectin 9, trim16, and trim25) in anterior kidney of salmon by lice infection. In addition, these authors found increased infectious salmon anemia viral load and virus-mediated mortality in L. salmonis-infected salmon [13]. These findings imply that antiviral response was compromised during the lice-infection.

Miscellaneous Elements in Innate Immunity
A panel of genes involved in innate immunity was downregulated in ATT and ADJ compared to PRE (segment 1; hdd11, cd99, c1qtnf6) or upregulated in ATT compared to PRE (segment 2; fgl1, camp, arg2, nlrc3) (Table 2). HDD11, a protein of unknown function that was bacteria-induced in silkworm [100], showed remarkable transcriptional suppression (~25-fold) in the present study. Downregulation of cd99 indicates that the regulation of lymphocyte adhesion [101] might be affected in chalimus-attached fins. CTRP6, encoded by c1qtnf6, induces IL10 and is considered to be anti-inflammatory [102,103]. While these modulated transcriptional patterns indicated a favorable environment for parasitism, the expression of some other genes suggested that salmon mount certain defense measures. For instance, FGL1 is an immunosuppressive ligand of the LAG3 receptor that plays a T cell inhibitory role [104]. The mRNA level of fgl1 was significantly suppressed in our study (Table 2). Two camp gene copies were found in Atlantic salmon and rainbow trout [105]. In the present study, these two paralogs demonstrated distinct expression profiles. Based on microarray, both camp paralogs were upregulated in ATT compared to ADJ, however this was only confirmed by QPCR for camp-a (Figures 6 and 8A,E). Mitochondrial expression of ARG2 regulates nitric oxide synthesis, and its transcript upregulation, as reported earlier in salmonid skin [39], was suggestive of M1 (classically activated) macrophage infiltration into the cutaneous region at chalimus-attachment sites. NLRC3, that plays inhibitory roles during inflammation by impeding NF-κB activation [106], was significantly upregulated at the mRNA level in response to lice in ATT and ADJ compared to PRE ( Table 2). These results indicate that chalimus-attachment modulated various innate immune elements in salmon fin.

Adaptive Immunity
Our microarray data identified different dysregulated elements of adaptive immunity, as well (e.g., two MH class I transcripts (hlab and hlah), b2m; Table 2). Downregulation of MH class I and II antigen genes at the lice-attachment site in skin was previously evidenced in salmonids [24,38,39]. Infiltrated MH class II+ cells in dermis and epidermis of lice-infected coho salmon, which demonstrated marked upregulation of mh class II in the skin, provided strong evidence for the contribution of MH -associated T cell-mediated immunity to the lice-resistance [39]. Additionally, the gene encoding B2M, a structural component of MHC I complex, was also found to be less abundant in ATT and ADJ compared to PRE, as observed in an earlier study [23].
Similarly, our data indicated decreased transcript levels of some B cell markers (iglc3, fcer1g) in infected fish. Different host-parasite models have earlier suggested that the immunosuppression caused by parasites is at least in part due to the decreased adaptive immune responses. During amoebic gill disease (AGD), MH class I pathway-related genes were downregulated in S. salar [107]. Suppression of mh class II mRNA was also observed in carp infected with Trypanoplasma borreli [108]. During the bi-phasic defense response in Atlantic salmon against sea lice, elevated levels of igm and igt transcripts in skin occurred only at mid-(2 weeks post-infection) and later-stages, suggesting impaired B cell-mediated immunity during the early stage of infection [24].

Relationship Between Transcript Expression and Lice Load
A clear separation between control and lice-infected groups was observed in the PCoA and hierarchical clustering analyses performed using either microarray or QPCR data. In contrast, ADJ and ATT groups showed overlap in the PCoA and a mixed distribution in clustering analyses ( Figure 1C,D and Figure 11). Our results suggest that contemporaneous regulation of elements responsible for opposing immune regulation was a part of the anti-lice response. For instance, several transcripts encoding positive regulators of immunity (e.g., cytokines: ccl17, ifna-a; APP: crp; redox markers: txn-b, prxl2a, prdx1a) demonstrated a negative correlation with lice burden ( Figure 11B). In contrast, transcript abundance of cd83, a well-known marker for mature dendritic cells and an immune-suppressive molecule [109], showed a significant positive correlation with lice load. These observations confirmed that lice-induced immune suppression in fin is one of the strategies that could negatively impact the susceptible hosts.

Compromised Immune System by Lice
In resistant species (coho and pink salmon), skin responses against sea lice infection are characterized by filamentous cell proliferation, moderately increased intracellular space, abundant leukocyte infiltration, epidermal thickening, and hyperplasia [36,64]. The magnitude of both inflammatory and hyperplasic responses was postulated as the primary determinant of resistant versus susceptible traits against sea lice infection in salmonids [10]. In addition, it has been found that excretory and secretory products of sea lice that contains immunomodulatory substances (e.g., trypsins and prostaglandins) could also have regulatory potential on host immunity [110]. Data presented in our current study provide evidence for hypo-inflammatory responses, weak APR, and compromised immunity by transcriptomic suppression of several components of host defense, such as antiviral responses, wound repair mechanisms, immune signaling pathways, and redox homeostasis during lice infection.
In an experimental infection trial, all the parasites will likely be in similar/closer developmental stage(s) of their life cycle and potentially demonstrate a development stage-specific impact on their host. However, in the natural environment, the infection dynamics will be more complex, with multiple life-stages of sea lice potentially parasitizing the host. The host-response would also likely be more complicated and depend on multiple factors. With this in mind, chronic lesions and open wounds resulting from lice infection and feeding could act as 'ports of entry' to potential secondary pathogens. In addition, the potential of L. salmonis as mechanical vectors in transmitting bacterial and viral diseases has already been experimentally demonstrated [111,112]. This phenomenon adds an extra layer of complexity to the host-pathogen interactions and opens a new paradigm, namely co-infection, which has been insufficiently studied until recently [13]. A compromised host immune system resulting from sea lice infection could be the main cause of detrimental consequences of co-infection [13]. Hence, exploring potential control measures (e.g., modulating host immunity using novel feed formulations [14]) to tackle sea lice infections, and potential secondary infectious outbreaks, requires further research attention. Molecular biomarkers identified herein will be valuable tools in our future endeavors.

Ethics Statement
All procedures involving fish handling, treatment, euthanasia, and dissection were performed in accordance with the guidelines of the Canadian Council of Animal Care (approved Memorial University Institutional Animal Care Protocol 17-77-MR, May 2, 2017).

Experimental Animals
Salmon smolts were purchased from a regional farm (Stephenville, NL, Canada) and transported to the Dr. Joe Brown Aquatic Research Building (JBARB, Ocean Sciences Centre (OSC), Memorial University of Newfoundland (MUN), Canada), where they were transferred into 3800 L tanks. After arrival, salmon were intraperitoneally PIT (passive integrated transponder)-tagged. Eighty post-smolts (307.3 ± 27.8 g mean initial weight ± SE) were transferred to the bio-containment zone in the Cold-Ocean Deep-Sea Research Facility (CDRF, OSC, MUN) and randomly distributed to and maintained under a 24 h light photoperiod in two 620 L tanks (40 fish/tank and 27 ± 1.5 kg/m 3 ) with a flow-through seawater supply. Fish were fed nightly with a commercial diet (Dynamic S, EWOS, Cargill) to satiation using automatic feeders. Fish were acclimated to laboratory conditions for four weeks before sea lice infection. Water quality parameters, such as temperature and oxygen saturation, were stably maintained throughout the experimental period (at 10 ± 1 • C and >90% O 2 saturation; Supplementary Tables S1 and S2). Fish were starved for 24 h before any handling or sampling.

Sea Lice Infection and Sampling
Adult female L. salmonis sea lice were collected from marine aquaculture sites in the Bay of Fundy (St. Andrews, NB, Canada) and transported to the Huntsman Marine Science Centre (St. Andrews, NB, Canada). Egg-strings were removed from adult females and maintained in hatching chambers at 10-12 • C for 9-10 days until they developed to infective free-living copepodid stage lice, which were then transported to CDRF.
Just before the infection, the water supply for tanks was turned off, the water level was reduced by 50%, and water was oxygenated with air diffusers to maintain~100% O 2 saturation at 10 ± 1 • C. Atlantic salmon were then bath-infected with L. salmonis by adding~50 copepodids per fish to each tank. Water quality parameters (dissolved O 2 and temperature) were continuously monitored throughout the infection period at intervals of 10 min. Water supply was reinstated after 2 h of copepodid exposure. At least five individuals from each tank at 0 days post-infection (dpi; prior to infection, PRE), 8 dpi, and 30 dpi (when most lice were at their chalimus and pre-adult stages, respectively) were euthanized with a lethal dose of MS-222 (400 mg/L; Syndel Laboratories, BC, Canada) and dissected for sample collection ( Figure 12A). Parasites on each salmon were counted. Samples from pelvic fin tips at an area of the fin with lice 'attached' (ATT) and 'adjacent' to chalimus attachment sites with no signs of previous lice attachment (ADJ) were excised as shown in Figure 12B, chalimi were removed from ATT sites using sterile forceps, and fin samples were flash-frozen in liquid nitrogen and stored at −80 • C until RNA extraction was performed.

RNA Extraction, DNA Digestion, and Column Purification
Fin tips were homogenized in TRIzol ® Reagent (Invitrogen/Life Technologies, Burlington, Canada) using a TissueLyser II (QIAGEN, Restch, GmbH, Germany), and total RNA was extracted from the homogenate as per the manufacturer's instructions. To eliminate any residual genomic DNA contamination, 40 µg of crude RNA from each sample was treated with DNase I (6.8 Kunitz units; RNase-Free DNase Set, QIAGEN) and then column-purified by using the RNeasy Mini Kit (QIAGEN) following the manufacturer's protocol. The RNA concentration and purity were assessed by ND-1000 UV spectrophotometry (NanoDrop, Wilmington, DE, USA), and the RNA integrity was examined by 1% agarose gel electrophoresis. The samples with tight 18S and 28S ribosomal RNA bands and high A260/280 and A260/230 ratios (>1.8) were used in transcriptional analyses.

Microarray Experiment: Design, Hybridization, and Data Acquisition
The microarray experiment was designed and conducted in compliance with the MIAME (Minimum Information About a Microarray Experiment) guidelines [113]. Fin RNA samples of six fish at 8 dpi (three from each duplicate tank) each from control (PRE) and sea lice-infected salmon (both ATT and ADJ samples collected from the same fish), selected based on RNA quality, were used in the microarray experiment (i.e., 12 animals and 18 samples in total). The microarray study was based on a common reference design ( Figure 12C), in which the differences among sample groups in fin transcriptome were investigated by contrasting individual samples against a common reference pool prepared from the equal contribution of all 18 samples. A consortium for Genomic Research on All Salmonids Project (cGRASP)-designed 44K Atlantic salmon oligonucleotide array (GEO accession: GPL11299; Agilent Technologies, Mississauga, Canada) that comprises 60mer probes was used as the microarray platform [27]. Pelvic fin tips were sampled from chalimus-attachment (ATT) and adjacent (ADJ) sites. (C) Common reference-based microarray experimental design. Arrows represent the biological replicates used from three treatment groups (PRE, fins from control fish prior to lice exposure (n = 6); ATT and ADJ, fins from chalimus-attachment and adjacent sites, respectively, both sampled from same fish (n = 6)). The base and arrowhead show the Cy3-labelled reference pool and Cy5-labelled experimental samples, respectively.
The antisense amplified RNA (aRNA) preparation and labeling were conducted using the Amino Allyl MessageAmp II aRNA amplification kit (Invitrogen). One microgram of DNase I-treated, column-purified total RNA was reverse-transcribed to the first-strand cDNA and used in second-strand cDNA synthesis. The double-stranded cDNA was then used to synthesize modified aRNA using the manufacturer's instructions. The aRNA samples were column-purified and subjected to the quantity and quality assessment by NanoDrop spectrophotometry and agarose gel electrophoresis, respectively. A common reference was prepared by pooling 10 µg of each of the 18 experimental samples. Twenty micrograms of each experimental and common reference aRNA sample were precipitated overnight through standard molecular biology methods and then resuspended in coupling buffer. While 18 individual experimental aRNA samples were labeled with Cy5, the common reference was labeled with Cy3 (GE Healthcare, Mississauga, ON, Canada) through a dye-coupling reaction, as per the manufacturer's instructions. The labeling efficiency was determined using the 'microarray' mode in the ND-1000 NanoDrop spectrophotometry. For each array, equal amounts (825 ng) of each Cy5-labelled experimental aRNA and Cy3-labelled reference aRNA were pooled, fragmented following the manufacturer's protocol (Gene Expression Hybridization Kit; Agilent, Santa Clara, CA, USA), and co-hybridized to a 44K salmonid oligonucleotide microarray. Each array was hybridized at 65 • C for 17 h at a 10 rpm rotation in an Agilent hybridization oven. Array slides were then washed with Gene Expression Wash Buffer 1 and 2 (Agilent) according to the manufacturer's guidelines, and residual wash buffer was removed by centrifuging at 200× g for 5 min at room temperature.
Each microarray was scanned at 5 µm resolution using a SureScan D Microarray Scanner (G2600D, Agilent Technologies) using Agilent Scan Control Software (v9.1.11.7, Agilent Technologies) by applying a built-in protocol (Agilent_HD_GX_2color). Photomultiplier tube (PMT) sensitivity for Cy3 and Cy5 dye channels was set to 100%. The resulting TIFF images (20 bit) containing raw array data were then subjected to Agilent Feature Extraction (FE, v12.0.3.1) software to retrieve the signal intensity data of each microarray probe using a protocol provided by the developer (GE2_1200_June14). Probe signal, background signal subtraction, LOWESS normalization, outlier flagging, and log 2 -ratio calculations were performed by FE. Quality control (QC) parameters from QC reports of each array generated by FE were also examined. Resulting text files containing raw microarray data were then subjected to GeneSpring GX, a multi-omic analysis platform (v14.9.1, Agilent Technologies), to visualize and further process the data. The missing values were imputed using a method embedded within GeneSpring GX. Based on QC assessment and downstream analyses (such as heatmaps and principal component analyses (PCA)), samples considered as potential outliers were identified (1 fish/group) and excluded from the downstream analyses. Features absent in more than 25% of the arrays in all three groups (PRE, ADJ, and ATT) were discarded, resulting in a final list that consisted of 28,065 probes for statistical analyses. Microarray data from this study have been submitted to Gene Expression Omnibus (GEO) under the accession GSE140756.

Microarray Data Analysis
The DEPs between sample groups (PRE, ADJ, and ATT) were determined by Significance Analysis of Microarrays (SAM) [31] using the Bioconductor package siggenes in the R environment (v1.1.463). Two-class comparison with a false discovery rate (FDR) cutoff of 0.05 or 0.01 was performed to identify the chalimus-responsive probes between groups (i.e., ATT versus PRE, ADJ versus PRE, and ATT versus ADJ). The resulting transcript list identifications were updated using the contiguous sequences (contigs) or expressed sequence tags (ESTs) that were used to design the 60mer probes of the array [27]. BLASTx searches of these transcript sequences against the NCBI non-redundant amino acid (nr) and Swiss-Prot databases were performed with default settings (E-value threshold < 1 × 10 −5 ). The resulting BLASTx hits were mapped to gene ontology (GO) terms of chalimus-responsive transcripts in each sample group (biological process (BP), molecular function (MF), cellular component (CC)) using the Blast2GO program (BioBam Bioinformatics, Valencia, Spain) [114].
The annotation assessment of this salmonid 44K 60mer oligo microarray resulted in >10K probes with "unknown" annotations. In order to update the annotation, 60mer probes and their representative contig sequences were searched against the Swiss-Prot database (April 2019 version) and the NCBI nr/nt database. The resulting hits were filtered with E-value < 1 × 10 −5 , identity percentage (>75%) and query coverage percentage (>50%). The filtering criteria for 60mer probe BLASTn hits were stringent (only 2 allowed mismatches with un-gapped alignment option). The probe annotations were also revised by homology sequence searches in updated genome annotation databases for S. salar, O. mykiss, Danio rerio, and Homo sapiens. The gene symbols for probes were assigned from HUGO Gene Nomenclature Committee (HGNC; https://www.genenames.org/) and/or GeneCards (https://www.genecards.org/) databases.

GO Term Enrichment Test, Analyses and Visualization of GO/Pathway Term Networks and Expression Profiles
Initially, a GO term enrichment test was performed on DEP lists, which were subjected to a fold-change (FC) cutoff of |2|, by using the Blast2GO (Fisher's Exact Test, FDR cutoff of 0.01) with the 44K salmon array as the reference set. Secondly, the functional implications of the lists of differentially expressed genes (DEGs) were examined using ClueGO [32] plugin in Cytoscape (v3.5.1) [33]. The GO databases (23.05.2019) for BP, MF, and CC, and the Reactome pathways database (23.05.2019) were used. The enrichment/depletion analysis was performed using a two-sided hypergeometric test after its adjustment by the Bonferroni step-down procedure. The GO term fusion strategy was employed to integrate GO categories, minimize the complexity, and create a functionally organized GO cluster network. The leading terms were ranked based on their significance with different p-values for different gene sets with Kappa-statistics score threshold set to 0.4.
We performed a Pathifier pathway analysis that incorporated transcriptional expression data to infer dysregulated pathways. All the significant GO/pathway terms from the above section were pooled and their Pathway Deregulation Score (PDS) based on expression was computed by using Pathifier Bioconductor package 1.16.0 in the R environment. This method by Drier et al. [34] calculates the PDS for each pathway per sample by using the annotation resources such as the MSigDB (http://software.broadinstitute.org/gsea/msigdb. Using the C5 collection of MSigDB, 134 GO/pathway terms were annotated and the PDS was estimated for infected salmon (ADJ and ATT) with respect to the non-infected control (PRE).
Pearson correlation and complete linkage clustering functions in the Genesis program were used in generating the hierarchical clustering and heatmaps of median-centered data of DEPs [115]. Alternatively, unsupervised clustering on Euclidean distance of PDS values was constructed using Ward's agglomerative linkage method (ward.D2) and visualized as a heatmap by hierarchical clustering using the heatmap3 R package [116].

Real-Time Quantitative Polymerase Chain Reaction (QPCR) Confirmation
Thirty-one microarray-identified GOIs, including some paralogs for 5 GOIs, were selected from the gene lists for transcriptional profiling by QPCR. GOIs were selected based on the functional themes they represented (i.e., 'ECM catabolism, tissue remodeling, and wound healing', 'immunity and defense', 'chemotaxis and signaling', 'antiviral response', and 'redox homeostasis'; Supplementary  Table S3). Another member of the MMP family (mmp20), which was not identified by the microarray analyses, was also added to the QPCR analysis.
In addition to the samples used in the microarray experiment, we included 4 more samples from each group (PRE (10 fish), ADJ and ATT (10 fish)) in the QPCR-confirmation experiment (n = 10/group; 30 samples in total). First-strand cDNA templates for QPCR were synthesized in 20 µL reactions from 1 µg of DNaseI-treated, column-purified total RNA using random primers (250 ng) and M-MLV reverse transcriptase (200 U; Invitrogen/Life Technologies) with the manufacturer's first strand buffer (1× final concentration), dNTPs (0.5 mM final concentration), and DTT (10 mM final concentration) at 37 • C for 50 min. The resulting cDNA was diluted either 20 times (for QC; see below) or 40 times (for normalizer testing and expression study) with nuclease-free water (Invitrogen). QPCR assays were performed in technical duplicates (for QC) or triplicates (for normalizer testing and expression study) using Power SYBR Green I dye chemistry in 384-well format on a ViiA 7 Real-Time PCR system (Applied Biosystems, Thermo Fisher Scientific). Each QPCR plate included a no-template control (NTC). Each assay was performed in 13 µL reaction containing 1× Power SYBR Green PCR Master Mix (Applied Biosystems/Life Technologies), 50 nM of both the forward and reverse primers, and 4 µL of diluted cDNA (5 ng input total RNA for expression study and normalizer testing, and a serially-diluted RNA for amplification efficiency determination; see below). The QPCR profile consisted of 1 cycle of 50 • C for 2 min, 1 cycle of 95 • C for 10 min, and 40 cycles of 95 • C for 15 s and 60 • C for 1 min, with fluorescence detection at the end of each 60 • C step.
Details of each QPCR assay, such as primer sequences, amplification efficiency, and amplicon size, are presented in Supplementary Table S3. BLASTn searches against each GOI were performed using nr/nt and EST databases of NCBI in order to identify any existing paralogs. By aligning the sequences of identified putative paralogs in Vector NTI (Vector NTI Advance 11, Life Technologies), potential regions (i.e., with at least 2 bp difference) for designing paralog-specific primers were identified. Primers were designed using either Primer 3 (v4.1.0; http://bioinfo.ut.ee/primer3/) or PrimerQuest (Integrated DNA Technologies, Inc.; https://www.idtdna.com/Primerquest/Home/Index). Prior to performing QPCR assays, each primer pair underwent a series of QC procedures. In brief, a reference RNA pool was prepared using four randomly selected RNA samples from each sample group (12 in total). Amplification efficiencies [117] were determined from the standard curves generated using a 5-point 1:3 dilution series starting with cDNA corresponding to 10 ng of input total RNA. It was ensured that each primer pair demonstrated a single peak in the dissociation curve analysis and had an amplification efficiency > 90% with no amplification in NTCs.
After completing the primer QC and normalizer test, QPCR analyses were performed in compliance with MIQE (Minimum Information for Publication of Quantitative Real-Time PCR Experiments) guidelines [119] to determine the transcript expression of GOIs. The relative quantity (RQ) of each transcript was determined by a qBase relative quantification framework [120,121] by using the C T values measured for GOIs and reference genes by the ViiA 7 Software (v1.2.3; Applied Biosystems, Thermo Fisher Scientific). RQ of each GOI was calculated by normalizing against transcript levels of both ef1a1 and eif3d, and by incorporating the amplification efficiencies. The individual with the lowest normalized expression level was considered as the calibrator (i.e., RQ = 1) for each GOI.

Statistical Analyses
RQs of GOIs were checked for outliers using Grubb's test, and for homoscedasticity using Levene's test. In total, 21 RQ values were identified as statistical outliers in the entire dataset (i.e., out of 1110 RQ values), and were excluded from the study. The differences among treatments were then compared using one-way ANOVA followed by Tukey's post hoc test. The accepted level of significance was p < 0.05. All the statistical analyses above were conducted using IBM SPSS Statistics (v25.0.0; IBM Corp, Armonk, NY).
FC values derived from microarray and QPCR were log 2 -transformed and analyzed for correlation via linear regression. A significant correlation between both datasets was considered as evidence for the validity of the microarray results.
Normalized log 2 ratios from microarray data were transformed (by adding a constant so that values were positive) and analyzed using Principal Coordinates Analysis (PCoA) for similarities among the expression patterns of DEPs observed among the salmon. Similarly, standardized RQ values of GOIs from QPCR were used along with phenotypic parameters (i.e., chalimus count, weight, length, hepatosomatic index (HSI), spleen-somatic index (SSI), viscerosomatic index (VSI), and condition factor (K)) in PCoA. We also tested for differences among gene expression patterns via permutational multivariate analysis of variance (PERMANOVA) with 9999 random permutations. The accepted level of significance was p < 0.05. Both PERMANOVA and PCoA were based on Bray-Curtis similarities of all pairwise comparisons among individuals. These analyses were performed using PRIMER 6.1.15 (Ivybridge, UK). The complete linkage hierarchical clustering was conducted with PRIMER using Pearson correlation resemblance matrices. Log 2 -transformed phenotypic and QPCR-gene expression data were subjected to Pearson correlation analyses using IBM SPSS to find any relationship between them, and the resulting Pearson correlation coefficients (r) were tabulated as a correlation matrix.

Conclusions
Using a laboratory infection model and molecular techniques (i.e., microarray and QPCR), we investigated and profiled the transcriptome response in fin tissue at the early (chalimus) stage of sea lice infection in Atlantic salmon. Findings from the current study significantly contribute to the understanding of the physiological basis for Atlantic salmon's high vulnerability to sea lice infection, particularly at the chalimus stage. Immunosuppression at transcriptional levels of several GOIs in ATT and ADJ sites compared to PRE suggested lice-associated immunomodulation, by which chalimus disrupts or evades the defense strategies of the host. A comparison between ATT and ADJ sites, however, indicated that Atlantic salmon attempt to mount a local anti-parasite response at the site of chalimus-attachment. Transcriptional changes that were documented in this study in response to chalimus infection on fin tissue may help to explain the Atlantic salmon's inability to mount a tissue response against sea lice that is sufficiently robust to expel the parasite.