Integration of Transcriptome, Gross Morphology and Histopathology in the Gill of Sea Farmed Atlantic Salmon (Salmo salar): Lessons From Multi-Site Sampling

The gill of teleost fish is a multifunctional organ involved in many physiological processes such as gas exchange, osmotic and ionic regulation, acid-base balance and excretion of nitrogenous waste. Due to its extensive interface with the environment, the gill plays a key role as a primary mucosal defense tissue against pathogens, as manifested by the presence of the gill-associated lymphoid tissue (GIALT). In recent years, the prevalence of multifactorial gill pathologies has increased significantly, causing substantial losses in Atlantic salmon aquaculture. The transition from healthy to unhealthy gill phenotypes and the progression of multifactorial gill pathologies, such as proliferative gill disease (PGD), proliferative gill inflammation (PGI) and complex gill disorder (CGD), are commonly characterized by epithelial hyperplasia, lamellar fusion and inflammation. Routine monitoring for PGD relies on visual inspection and non-invasive scoring of the gill tissue (gross morphology), coupled with histopathological examination of gill sections. To explore the underlying molecular events that are associated with the progression of PGD, we sampled Atlantic salmon from three different marine production sites in Scotland and examined the gill tissue at three different levels of organization: gross morphology with the use of PGD scores (macroscopic examination), whole transcriptome (gene expression by RNA-seq) and histopathology (microscopic examination). Our results strongly suggested that the changes in PGD scores of the gill tissue were not associated with the changes in gene expression or histopathology. In contrast, integration of the gill RNA-seq data with the gill histopathology enabled us to identify common gene expression patterns associated with multifactorial gill disease, independently from the origin of samples. We demonstrated that the gene expression patterns associated with multifactorial gill disease were dominated by two processes: a range of immune responses driven by pro-inflammatory cytokines and the events associated with tissue damage and repair, driven by caspases and angiogenin.


INTRODUCTION
The gill of teleost fish is a multifunctional organ involved in many physiological processes such as gas exchange, osmotic and ionic regulation, acid-base balance and excretion of nitrogenous waste (Evans et al., 2005). To facilitate these functions, the gill tissue has evolved into a highly complex system of branching vascular structures (lamellae), separated from the external milieu only by a thin layer of gill epithelium and mucosa (Koppang et al., 2015;Salinas, 2015). The densely packed lamellar structure of the gill is highly advantageous because it provides a large surface area for oxygen transfer, amounting to approximately 0.1-0.4 m 2 of lamellar surface per kg of body mass (Maina, 2011;Park et al., 2014). However, having such an extensive interface with the environment comes at a high price (reviewed in Nilsson, 2007). Firstly, the large respiratory surface area may contribute to the increased water and ion fluxes that need to be counteracted by energetically expensive ion pumping. That so-called osmorespiratory compromise, which is a trade-off between high gill permeability (to promote respiratory gas exchange) and low gill permeability (to limit unfavorable water and ion fluxes), has been demonstrated in freshwater fish during exercise and in hypoxia-intolerant species (e.g., salmonids) during the periods of exposure to low-oxygen water (Onukwufor and Wood, 2018;Morgenroth et al., 2019). Secondly, the highly complex vasculature of the gill is prone to mechanical injuries, increasing the risk of hemorrhage. Because the gill is the only organ in fish that receives the entire cardiac output, large or frequent gill hemorrhages can be fatal. Thirdly, the large lamellar surface area may facilitate the uptake of toxic substances, both those that occur naturally (such as ammonia, algal toxins and metal ions) as well as man-made pollutants (e.g., pesticides, herbicides, fertilizers, detergents, industrial chemicals and microplastics) (McIntyre et al., 2018;Kintner and Brierley, 2019). Finally, the close proximity of blood and environmental water across the large respiratory surface area provides (1) major ports of entry for pathogens via transepithelial transport and (2) attachment sites for water-born parasites. To combat pathogens and parasites, the gill is equipped with the gill-associated lymphoid tissue (GIALT), adding another dimension to the importance of this organ to fish functioning and survival (Koppang et al., 2015;Salinas, 2015;Xu et al., 2020).
To avoid 'the large interface' penalties, fish constantly adjust the functional size of the gill tissue to match their current oxygen demands (Nilsson, 2007). When faced with changes in their need for oxygen uptake (due to changes in metabolic rate or water oxygen content), fish are known to (1) regulate the water flow over the gill by adjusting the volume and frequency of buccal pumping and/or (2) regulate the blood flow inside the gill to alter the perfusion levels of lamellae. New evidence suggests that many teleost species may also undergo extensive gill tissue remodeling, during which lamellae are either completely embedded in a cell mass (to minimize the respiratory surface area) or fully protruded (to maximize the functional size of the tissue), with many potential intermediate stages between these extremes (Sollid et al., 2003;Anttila et al., 2015). Whether these regulatory and compensatory mechanisms aiming to optimize gill performance are sufficient for fish to cope with climate change, especially in the areas of high human impacts, remains unknown.
Fish are predominantly ectothermic and any increases in water temperature resulting from climate change will increase their metabolic rate and thus demands for oxygen. The Q 10 for metabolism in fish (the fold increase in metabolic rate with a 10% increase in water temperature) varies from 1.5 to 2 (White et al., 2006). But the elevated oxygen demand is not the only challenge that the gill tissue is facing. Warmer water holds less oxygen and rising global temperatures alter oceanic circulation, which altogether may decrease the oxygen supply available for fish to support their increased metabolism (McBryan et al., 2013). It has been estimated that the open ocean lost ∼2% (77 billion metric tons) of its oxygen over the past 50 years (reviewed in Breitburg et al., 2018). The extent of deoxygenation of coastal waters may even be greater, due to eutrophication from agriculture and sewage runoff (Diaz and Rosenberg, 2008). Furthermore, environmental changes, including temperature increases, have been linked to enhanced expression of marine infectious diseases (Burge et al., 2014). Many marine organisms, including marine pathogens, are shifting their distributions poleward or to deeper waters as ocean temperatures warm (Nye et al., 2009;Pinsky et al., 2013). As a result, many species of fish are at risk of being exposed to novel pathogens and parasites that were not part of their evolutionary history (LeBlanc et al., 2019). Climate change may also affect the fish-pathogen interactions that are already well established, for example by altering pathogen virulence, overall pathogenicity and/or host immune response (Benedicenti et al., 2019;Laurin et al., 2019). Given the fact that in fish most environmental stressors converge at the level of the gill, performance of this tissue and its capacity to accommodate environmental change are central to understand short-and long-term impacts of warming and hypoxia on fish species (Akbarzadeh et al., 2018;McCormick and Regish, 2018;Houde et al., 2019). Monitoring gill health is especially important in farmed marine fish, which are restricted to coastal waters and have limited microhabitat choice.
In recent years, the prevalence of gill damage and disease in farmed marine fish has significantly increased, leading to substantial losses in the Scottish and global Atlantic salmon (Salmo salar) aquaculture industry (Young et al., 2007;Shinn et al., 2015;Nowak and Archibald, 2018). Although gill-related mortality varies from site to site and between years, July to December in Northern Hemisphere is considered the higher risk period because the sea during these months is warmer. With the projected increase in sea water temperature associated with climate change, gill pathologies are now regarded as the greatest challenge to the farmed salmon sector due to direct and indirect impacts on productivity. Losses may occur through direct mortalities, poor growth rates linked to decreased feed efficiency and because of the increased risk of co-infection (Rodger, 2007;Downes et al., 2018). The transition from healthy to unhealthy gill tissue happens through the progression of various gill pathologies such as amoebic gill disease (AGD), proliferative gill disease (PGD), proliferative gill inflammation (PGI) and complex gill disorder (CGD), which are commonly characterized by epithelial hyperplasia, lamellar fusion and inflammation (Steinum et al., 2010;Matthews et al., 2013;Mitchell et al., 2013;Herrero et al., 2018;Gjessing et al., 2019). The diagnostic of gill pathologies is often difficult because single-cause diseases (such as AGD) are compounded by multi-cause (multifactorial) disorders such as PGD, PGI or CGD, with complex etiology and potentially synergistic effects of co-infections (Gjessing et al., 2017;Gunnarsson et al., 2017;English et al., 2019). Due to overlapping symptoms of multifactorial gill inflammation, the terms PGD, PGI and CGD are often used interchangeably, with the first two historically referring to salmon farmed in Scotland (PGD) and Norway (PGI), while the last one highlighting the growing complexity of gill diseases and encompassing the syndromes referred to as PGD and PGI (Rozas-Serri, 2019). The multifactorial gill diseases are associated with a number of causative agents, including phytoplankton, zooplankton, viruses, fungi, bacteria and larger parasites (Gjessing et al., 2017;Herrero et al., 2018), with climate change contributing to emerging new aquaculture pathogens (Bayliss et al., 2017). Furthermore, gill diseases not only directly affect fish health and performance, but also creates a challenge for farmed salmon producers as they are required to treat compromised fish, increasing the risk of mechanical injuries to the gill vasculature (Powell et al., 2005). To meet the challenge of salmon gill health, the sector needs non-invasive monitoring tools for (1) early detection of gill pathologies, (2) diagnostics that is sufficient to establish appropriate treatment and (3) evaluation of effectiveness of treatment and recovery of gill function. It is important that such initiatives target multifactorial diseases, as they are becoming progressively more prevalent (Rahel and Olden, 2008).
The current practice is to regularly check sea farmed salmon for any signs of gill inflammation by in situ gross morphology examination of gill arches, using the PGD scoring system from 0 (no macroscopic pathology) to 5 (severe macroscopic pathology) (Gill Health Initiative, 2015;Bloecher et al., 2018). This is a non-invasive procedure that requires only light anesthesia and can be performed on a weekly basis in fish from across a facility. However, the applicability of the PGD scoring system to grade inflammation has never been tested in a systematic way. Because PGD is a multifactorial disease, progression from low to high PGD scores is likely to have a large component of case specificity, depending on the combination of pathogens involved, the overall health status of the fish and the environmental conditions. Despite complex etiology, the phenotypes generated by PGD are macroscopically similar (inflammation and hyperplasia of respiratory epithelium), which suggests that at least some of the underlying mechanisms associated with this pathology may be common (Mitchell and Rodger, 2011;Bruno et al., 2013). Unlocking these shared mechanisms is essential for early detection of gill disease, developing strategies to improve gill health and for fine-tuning salmon husbandry practices to the challenging conditions of a rapidly changing marine environment.
To identify the common gene expression patterns associated with gill inflammation, we examined Atlantic salmon from three marine production sites in Scotland. The sites and time of sampling were chosen to ensure (1) differences in fish cohort (hatchery origin, on-site treatments and overall health status), (2) diversity of pathogen and (3) differences in local environment. For each site, we selected fish with low and high PGD scores, analyzed their gill transcriptome (RNA-seq) and gill histopathology (microscopic examination), and then integrated these data to explore gene expression patterns associated with multifactorial gill disease.

Sampling and Gross Morphology
Fish (Atlantic salmon, Salmo salar) were sampled at three marine production sites belonging to Scottish Sea Farms (A on Isle of Mull and B and C in Shetland) between October 2017 and March 2018 (for details see Table 1 and Supplementary Figure 1). All fish were of strain Fanad and originated from the same egg fertilization batch. They were reared in different hatcheries (Couldoran, Pettigo-Damph and Knock-Frisa for sites A, B and C, respectively) for one year and entered the sea in spring 2017.
During sampling, randomly selected fish (76 in total) were placed in an anesthetic bath (∼20 g of MS-222/150 L) for 5-10 min, followed by macroscopic examination (gross morphology) of gill tissue as per routine procedure to monitor fish welfare indicators. This procedure is performed by fish health professionals, using semi-quantitative 6-grade scoring systems for PGD, from 0 (no visual pathology) to 5 (severe visual pathology) (Gill Health Initiative, 2015;Bloecher et al., 2018). Each of the 8 gill arches (4 on left and 4 on right side) was scored separately, generating in total 8 PGD scores per fish. Although scoring for PGD is typically non-destructive and requires only light anesthesia (e.g., 10 g of MS-222/150 L), fish in our study were subjected to terminal anesthesia for subsequent harvesting of tissue samples. Immediately after PGD scoring, fish were bled and the gill arch with the highest PGD score was excised for both transcriptome profiling and histopathological examination. For transcriptomic profiling, 3 transverse sections (gill arch ∼3 mm long with 3-4 pairs of filaments) from dorsal, medial and ventral regions of the gill (Figure 1) were transferred to RNAlater R (Sigma-Aldrich, St. Louis, MO, United States), kept at  (Cadoret et al., 2013) for 24 h, followed by a short-term exposure to 10% neutral buffered formalin before tissue processing for histopathological examination.

Selection of Fish With Low and High PGD Scores
Grouping all 76 sampled fish by their highest PGD score and site demonstrated that each site had fish with PGD scores 1, 2 or 3, but not 0 or 5, and only rarely 4 ( Table 2). Thus, the comparison between fish with low and high PGD scores was restricted to PGD scores 1 and 3, respectively. In total, 45 gill arches (27 with PGD score 1 and 18 with PGD score 3) were subjected to total FIGURE 1 | Approach to representatively divide the gill arch with the highest PGD score between transcriptome profiling and histopathological examination. For transcriptomic profiling, 3 transverse sections from dorsal, medial and ventral regions of the gill were subjected to total RNA extraction to generate individual RNA samples, which were subsequently pooled. All remaining tissue was used for histopathology.
RNA extraction, but the final RNA-seq sample size dropped to 43 gill arches (from 43 fish), due to losses at the both pre-and post-sequencing stage ( Table 2).

RNA Extraction
Total RNA extraction was performed on individual gill transverse sections (45 gills × 3 sections) after removing the arch tissue and leaving only full-length filaments for further processing (Figure 1). Briefly, the RNA was isolated by homogenization of ∼100 mg of gill tissue in TRIzol R Reagent (Ambion by Life Technologies, Carlsbad, CA, United States), using 3 mm tungsten carbide beads and a TissueLyser II Disruption System (Qiagen GmbH, Hilden, Germany). Following isolation, the RNA was quantified by spectrophotometry (NanoDrop Technologies, Wilmington, DE, United States) and its integrity was confirmed by electrophoresis (Agilent Technologies, Santa Clara, CA, United States). The three individual RNA samples that originated from the same gill were then pooled to generate a single RNA sample per gill tissue (n = 45 RNA pools in total), with an equimolar contribution of RNA from dorsal, medial and ventral regions of the gill to each pool. All but one pooled gill RNA samples had a 260/280 ratio > 1.8 and RIN number > 9.3, thus meeting the criteria for RNA-sequencing. The sample with degraded RNA was eliminated from further processing.
multi-mapping read, with n reflecting the total number of alignments reported for that read.

Differential Gene Expression
Differential expression analysis was performed using the Bioconductor package edgeR v3.22.5 (Robinson et al., 2010) in R v3.5.1 (R Core Team, 2018). Genes with a CPM (count per million) < 1 in three or more samples were removed, resulting in 35996 genes for analysis. Filtered counts were subsequently normalized using a trimmed mean of M-values (TMM) between each pair of samples. Based on exploratory data analysis, one library was identified as an outlier and removed from subsequent analysis (for details on sample size see Table 2). Data for the remaining 43 fish were modeled using a negative binomial generalized log-linear mixed model that included both group (PGD scores 1 and 3) and site (A, B and C) as fixed effects. In total, four contrasts were generated from the same model: PGD 3 vs PGD 1 (with the site effect blocked) and then sites A vs C, B vs C and A vs B (with the group effect blocked). Differentially expressed genes (DEGs) were identified at false discovery rate (FDR) < 0.01 and absolute Log 2 fold change (FC) > 1.

Functional Analysis of Gene Expression
Salmon DEGs were mapped to human orthologs to generate HGNC (HUGO Gene Nomenclature Committee) gene identifiers for functional analysis of the RNA-seq results. This approach has been demonstrated to improve biological interpretation of the salmon gene expression profiles by providing access to well-annotated databases and tools for mammalian model organisms, despite limitations of the mapping due to the extra genome duplication events in teleost fish and species-specific differences in gene function and molecular pathways (Song et al., 2014;Król et al., 2016). Mapping was done by aligning salmon transcript sequences to the protein sequences from the human genome (release 88, downloaded from Ensembl at 2 ) using BLASTX (version 2.2.31) (Camacho et al., 2009) with an E-value cut off of 0.00001 and a maximum of one target sequence for every transcript. As one transcript can have multiple hits against one target sequence, custom Python scripts were used to filter the blast results to contain only the hit with the highest identity for each transcript. Although the majority of the salmon genes mapped to a unique human ortholog, some salmon genes mapped to the same human ortholog (for details see Results). To obtain a single expression value (mean Log 2 FC) per HGNC gene identifier, the expression of the salmon genes mapped to the same human ortholog was either averaged (if contributing salmon genes had similar expression profiles) or based on the expression of the salmon gene that was more abundant (if contributing salmon genes had contrasting expression profiles). For biological interpretation of the RNA-seq results, human orthologs of the salmon DEGs along with their Log 2 FC values were analyzed using Ingenuity Pathway Analysis (IPA, QIAGEN Redwood City, www.qiagen.com/ingenuity) to explore (1) enrichment of canonical pathways, (2) upstream regulators and (3) downstream effects associated with these genes. The same set of genes was also submitted to PANTHER Classification System (Mi et al., 2013) to perform Gene Ontology (GO) enrichment analysis.

Gill Histopathology
Gill tissue was routinely dehydrated in ethanol, equilibrated in xylene and embedded in paraffin wax according to standard histological techniques (Bancroft and Gamble, 2007). Sagittal sections (3 µm) of the gill arch were cut with a microtome and mounted onto microscope slides. These sections were then subjected to haematoxylin and eosin (H&E) staining. All sections were digitized at 40 × magnification, using the Olympus dotSlide 2.1 Virtual Slide System (Olympus Corporation, Tokyo, Japan). The resultant images were randomized to ensure blinded examination and then scored using a system developed to assess gill histopathology in sea farmed Atlantic salmon (Mitchell et al., 2012), with minor amendments. The details of the scoring system are presented in Table 3.
Prior to the analysis, three variables were excluded (bE, bT and ob) as scores were invariant across fish (see Results). To compare the gill histopathology scores between PGD groups (1 and 3) and also between fish from different sites (A, B and C), non-metric multidimensional scaling (NMDS) was performed using the metaMDS function from the vegan package (Oksanen et al., 2019) in R. For the NMDS analysis, a Jaccards dissimilarity index was used to calculate the dissimilarity matrix, two dimensions were specified and 100 random starts used. Similarities between fish histopathology scores were visualized using a biplot with 95% confidence ellipses around the group centroids and variable vectors included. Other parasites or agents 1 Scoring system 0-3, with 0 = no pathological changes, 1 = mild changes affecting < 10% of gill tissue, 2 = moderate changes affecting 10-50% of gill tissue, and 3 = severe changes affecting > 50% of gill tissue.

Phenotypic Characteristics of Fish
Body mass of 43 fish subjected to the RNA-seq experiment varied from 1.5 to 4.0 kg, with their body length ranging from 47 to 62 cm (Figures 2A,B). The effects of group (PGD 1 and PGD 3) and site (A, B and C) on body mass and length were not significant (P > 0.05, 2-way ANOVA). As a result, the fish did not differ in their Fulton's condition factor K (P > 0.05, 2-way ANOVA), which varied from 1.1 to 1.7 ( Figure 2C). Given the small sample size ( Table 2), these results should be treated with caution.

Gill Transcriptome
The NMDS analysis of 43 gill transcriptomes revealed no differences between 26 fish with PGD score 1 and 17 fish with PGD score 3, as indicated by their overlapping 95% confidence intervals ( Figure 3A). This finding was reinforced by differential gene expression analysis, which showed 0 DEGs for comparison PGD 3 vs PGD 1 ( Table 4 and Supplementary Table 1). In contrast, gill transcriptomes were clearly separated by site, with no overlap in 95% confidence intervals for 20 fish from site A, 10 fish from site B and 13 fish from site C ( Figure 3B). The sitespecific differences in the gill transcriptomes are reflected in the number of DEGs identified between the sites: 1360 for A vs C, 708 for B vs C and 240 for A vs B (   Table 1.

Gill Histopathology
Histopathological examination of the gill tissue generated 15 scores (ranging from 0 to 3) per fish, focusing on 4 index and 11 ancillary criteria (Supplementary Table 5). Distribution of the histopathological scores by PGD score (1 and 3) and site (A, B and C) for each of 15 criteria is presented in Supplementary Figure 2. The NMDS analysis of the gill histopathology results for 43 fish showed no differences between 26 fish with PGD score 1 and 17 fish with PGD score 3, as indicated by their overlapping 95% confidence intervals ( Figure 4A). In contrast, gill histopathology differed between sites, with 1) fish from site C being clearly separated from other fish and 2) fish from site A overlapping with fish from site B ( Figure 4B). The position of the site-associated clusters along the NMDS1 axis indicates that fish from site C (with lower scores) show relatively low level of gill histopathology, while fish from sites A and B (with higher scores) have gill tissue with relatively moderate changes in histopathology. The main drivers of these differences are lamellar hyperplasia (LH), lamellar fusion (LF), cellular anomalies (CA) FIGURE 3 | Gill transcriptome of 43 fish grouped by PGD score (A) and site (B). Each panel shows a NMDS plot of gene expression profiles between different fish (numbers are fish IDs). The distances on the plots correspond to the leading fold change (FC), which is the average (root-mean-square) Log 2 FC for the 500 genes most divergent between each pair of fish. Ellipses indicate 95% confidence intervals, overlapping for fish grouped by PGD score (A) but not for fish grouped by site (B). The stress value of the NMDS ordination is 0.162.
and Neoparamoeba-like protist parasites (pp), as evidenced by the vectors paralleled to the NMDS1 axis in Figure 4B. Overall, fish from site A and independently from site B developed significantly higher grade inflammation than fish from site C, thus facilitating the search for the molecular bases of multifactorial gill disease.

Identification of Gene Expression Patterns Associated With Multifactorial Gill Disease
Our approach to identify the common gene expression patterns of non-specific gill inflammation is explained in Figure 5. Specifically, all gill transcriptomes from sites A and B (with moderate gill histopathology) were independently compared to all gill transcriptomes from site C (with low gill histopathology), which resulted in 1360 and 708 DEGs (at FDR < 0.01 and absolute Log 2 FC > 1) for A vs C and B vs C contrasts, respectively ( Table 4, for the lists of DEGs see Supplementary  Tables 2, 3). A relatively large number of genes were common between these two comparisons (462 DEGs in total), including 354 protein-coding genes, 35 immunoglobulin gene segments, 15 pseudogenes and 58 non-coding RNAs ( Table 5, for the list of common DEGs see Supplementary Table 6). The Log 2 FC values for these genes were averaged between A vs C and B vs C contrasts to provide one expression value per common gene. Importantly, all common DEGs show the same direction of change in both the moderate-to-low histopathology contrasts, i.e., they are either upregulated or downregulated in relation to site C, as indicated in Figure 5. That commonality of the response suggests that the identified DEGs are part of the gene expression profile indicative of progression from lower to higher grade inflammation of gill tissue. To understand the underlying mechanisms, we focused on the 354 protein-coding DEGs and predicted their functionality through mapping to the human orthologs.

Mapping Salmon Genes to Human Orthologs
The results of blasting all 462 salmon DEGs associated with multifactorial gill disease against the protein sequences from the human genome (BLASTX, E-value < 0.00001, top hit) are presented in Table 5 and Supplementary Table 6. Among them, 311 of 354 protein-coding salmon genes were matched to HGNC gene identifiers, but not all of them were unique. Specifically, 191 protein-coding salmon transcripts mapped uniquely to one human ortholog, while the remaining 120 protein-coding salmon transcripts mapped to 44 human orthologs, with > 1 salmon gene mapping to the same HGNC gene identifier (Supplementary Table 7). The expression of salmon genes mapped to the same  Table 3 and Supplementary Figure 2). Ellipses indicate 95% confidence intervals, overlapping for fish grouped by PGD score (A) and for fish from sites A and B but not from C (B). Vectors are fitted to visualize the contribution of histopathological parameters to the ordination. The stress value of the NMDS ordination 0.190.
human ortholog was averaged to provide one expression value (mean Log 2 FC) per HGNC gene identifier, if the response of the salmon genes was consistent (i.e., all contributing salmon genes were either upregulated or downregulated, but not both).
In two cases with the contrasting expression of salmon genes FIGURE 5 | Approach to identify gene expression patterns associated with multifactorial gill disease. Based on the histopathological examination, A and B were classified as two independent sites with relatively moderate gill histopathology, while C was classified as a site with relatively low gill histopathology. Firstly, the lists of differentially expressed genes (DEGs) for A vs C and B vs C contrasts were generated (FDR < 0.01 and absolute Log 2 FC > 1), yielding 1360 and 708 transcripts, respectively. Secondly, both lists of DEGs were checked for common genes (462 in total). Finally, the list of common genes was screened for a commonality of response to ensure that only genes either upregulated or downregulated in both sets of DEGs were considered to constitute a common gene expression profile of non-specific gill inflammation.  (Figure 6 and Supplementary Table 8).
Most of these pathways are associated with 1) cellular immune FIGURE 6 | Top canonical pathways altered in Atlantic salmon gills with moderate histopathology (sites A and B) vs low histopathology (site C). The analysis was performed on 311 salmon genes mapped to 235 human orthologs, using Ingenuity Pathway Analysis (IPA) and P-value < 0.01. For details and list of corresponding genes see Supplementary Table 8. response (IL-17 Signaling, IL-6 Signaling, Granzyme A Signaling, Crosstalk between Dendritic Cells and Natural Killer Cells, Agranulocyte Adhesion and Diapedesis and HMGB1 Signaling), 2) cytokine signaling (IL-17 Signaling, IL-6 Signaling, Acute Phase Response Signaling, Role of JAK family kinases in IL-6-type Cytokine Signaling, TNFR2 Signaling and HMGB1 Signaling) and (3) tissue damage and repair. The tissue damage and repair were evidenced by alterations of pathways related to cellular stress and injury (Autophagy and HMGB1 Signaling), apoptosis (TNFR2 Signaling and JAK/Stat Signaling), cellular growth (STAT3 Pathway and JAK/Stat Signaling) and proliferation and development (STAT3 Pathway and JAK/Stat Signaling).
IPA analysis of upstream regulators is based on prior knowledge of predictable effects between transcriptional regulators (e.g., transcription factors, cytokines, microRNAs, receptors, kinases, chemicals and drugs) and their target genes, stored in the Ingenuity Knowledge Base. When such analysis was performed on the gene expression profiles indicative of multifactorial gill disease, IPA identified 15 top upstream regulators associated with 101 of the 235 submitted DEGs, all highly significant (overlap P-value < 0.0001) and activated (z-score > 2) (Figure 7 and Supplementary Table 9). Among them were one endotoxin (LPS), one pattern recognition receptor (NOD2), 10 cytokines (IL1B, IFNG, TNF, IL2, LIF, IL12B, IL12A, IL1, IFNA2 and IL6), two growth factors (PDGF and AGT) and one transcription factor (NFKB).
IPA analysis of downstream effects predicts potential outcomes from gene expression data, using the Ingenuity Knowledge Base of differential gene expression in varying disease and functional states. The predictions made for gills with moderate histopathology (sites A and B) vs low histopathology (site C) included 13 downstream disease/functions at P-value < 9.4E-08, with nearly all DEGs (230 of 235) contributing to the effects (Figure 8 and Supplementary  Table 10 GO enrichment analysis identified 15 GO terms (Biological Process) that were overrepresented with the majority of DEGs (209 of 235) associated with multifactorial gill inflammation, generating fold enrichments from 1.2 to 2.6 at Bonferronicorrected P-value < 0.05 (Figure 9 and Supplementary  Table 11). The enriched GO terms pointed towards increased (1) intra-tissue communication (Cellular response to cytokine stimulus and Cell communication), (2) presence of external stimuli (Response to external stimulus and Response to oxygen-containing compound), (3) activated immune response (Immune system process) and 4) ongoing tissue remodeling (Animal organ development and Regulation of biological quality).

FIGURE 7 | Top upstream regulators predicted from gene expression in
Atlantic salmon gills with moderate histopathology (sites A and B) vs low histopathology (site C). The analysis was performed on 311 salmon genes mapped to 235 human orthologs, using Ingenuity Pathway Analysis (IPA). For details and list of corresponding genes see Supplementary Table 9.

Top Genes Associated With Multifactorial Gill Disease
Top genes were defined here as the DEGs with the largest magnitude of change in expression (absolute Log 2 FC > 2), thus including the protein-coding genes that were at least 4 times higher or 4 times lower expressed in gills with moderate histopathology (sites A and B) vs low histopathology (site C). The list of these top genes is presented in Table 6, with 43 upregulated salmon genes mapped to 25 unique human orthologs and 14 downregulated salmon genes mapped to 12 unique human orthologs.

DISCUSSION
Our study is the first to focus on the gill transcriptome of Atlantic salmon in a highly variable and complex environment of the marine production sites. We specifically chose to work on fish farmed in the open sea-based system (such as net pens), because (1) the use of net pens is currently the most common form of salmon production (Philis et al., 2019) and (2) in the open pens, only nets separate the fish from the environment, making them prone to developing gill pathologies but also assisting with the tissue recovery through high water exchange FIGURE 8 | Top downstream effects predicted from gene expression in Atlantic salmon gills with moderate histopathology (sites A and B) vs low histopathology (site C). The analysis was performed on 311 salmon genes mapped to 235 human orthologs, using Ingenuity Pathway Analysis (IPA). For details and list of corresponding genes see Supplementary Table 10. and oxygenation (Shinn et al., 2015;Nowak and Archibald, 2018). Farming in net pens is typically associated with the multifactorial gill diseases (PGD, PGI or CGD) rather than the single-cause gill pathologies, such as AGD (Herrero et al., 2018;Gjessing et al., 2019;Laurin et al., 2019). So far, the extensive transcriptomic profiling of the gill tissue has been performed only in salmon exposed to the controlled environment of the closed-system tanks, following a single-pathogen challenge to induce either AGD (Morrison et al., 2006;Wynne et al., 2008a;Bloecher et al., 2018;Boison et al., 2019;Robledo et al., 2019) or infectious salmon anemia (Valenzuela-Miranda et al., 2015). In most cases, fish were subjected to a single-dose challenge of the infectious organisms, with only few studies evaluating the effects of re-infection (Bloecher et al., 2018;Boison et al., 2019) or co-exposure to other infectious agents, such as hydroids (Bloecher et al., 2018) and Yersinia ruckeri (Valdenegro-Vega et al., 2015). To closely resemble the dynamics of multifactorial gill diseases observed in the production systems, fish need to be continuously exposed to a variety of infectious agents and environmental stressors that impact their health.
To evaluate the robustness and reliability of the PGD scores (gross morphology) in conveying information about gill health, we sampled salmon at three geographically distant locations in autumn and spring, without prior knowledge of the biotic and abiotic conditions of these sites. By doing this, we included in our data set the gill samples that were diverse in terms of origin, FIGURE 9 | Top GO terms (Biological Process) associated with differentially expressed genes (DEGs) in Atlantic salmon gills with moderate histopathology (sites A and B) vs low histopathology (site C). The analysis was performed on 311 salmon genes mapped to 235 human orthologs, using PANTHER overrepresentation test and Bonferroni-corrected P-value < 0.05. The GO terms are sorted by the fold enrichment of the most specific categories, with their parent terms indented directly below. For details and list of corresponding genes see Supplementary Table 11. time of sampling, fish background and the overall environmental milieu ( Table 1). The multi-site sampling approach to identify the common aspects of the multifactorial gill inflammation has been used before (Kvellestad et al., 2005;Steinum et al., 2010;Gjessing et al., 2019), but not in the context of the transcriptome analysis or underlying molecular mechanisms. Our study is the first to investigate the gill tissue from salmon at three different levels of organization: gross morphology (macroscopic examination), histopathology (microscopic examination) and whole transcriptome (gene expression) (Figure 1). We have focused on the gill arch with the highest PGD score per fish, as this arch may provide more accurate information about the general state of the gill health than the arch sampled randomly or the arch sampled because of its position, as practiced elsewhere (Wise et al., 2008;Steinum et al., 2010;Benedicenti et al., 2019;Gjessing et al., 2019). For better integration of the transcriptomic data with the gill gross morphology and histopathology, the tissue fragments used for the total RNA extraction aimed to represent the whole surface area of the gill rather than some specific areas of the tissue (Figure 1). This is in contrast to many gill gene expression (qPCR) studies, which tend to focus on either interbranchial lymphoid tissue or visible pathologies of the gill tissue, such as mucoid patches, hyperplastic lesions and lamellar fusions (Pennacchi et al., 2016;Marcos-López et al., 2018;Benedicenti et al., 2019).

PGD Scores as Measures of Gill Health
Effective treatment of salmon gill diseases requires development of diagnostic and prognostic tools that are non-invasive and suitable for a frequent use in the rapidly changing marine environment (Herrero et al., 2018;Gjessing et al., 2019;Rozas-Serri, 2019). We have specifically focused on the applicability of the PGD scores (gross morphology) to reflect gill health and underlying pathology, as this scoring system is already in use as part of routine fish welfare indicator assessment of salmon health in seawater farms worldwide (Gill Health Initiative, 2015;Bloecher et al., 2018). By comparing 43 gill samples with low (1) and high (3) PGD scores across three marine production sites in Scotland, we found that these two groups of gill tissue classified as different at the macroscopic level (PGD 1 and PGD 3) were in fact indistinguishable at the level of whole-transcriptome gene expression ( Figure 3A) and also indistinguishable at the level of microscopic histopathology ( Figure 4A). Furthermore, we could not identify any single gene that was expressed differently between the two groups ( Table 4).
Our results strongly suggest that the changes in gross morphology were not supported by the changes in gene expression or histopathology. The lack of detectable transcriptomic and/or histopathological changes associated with the progression of the PGD scores questions the suitability of the gross scoring system as diagnostic and prognostic tools to monitor and control both existing and emerging gill diseases. Instead, the PGD scores are good proxies for monitoring changes in macroscopic gill surface area available for respiration, the knowledge of which may be essential for fish health management and husbandry practices in aquaculture settings. The PGD scoring system of salmon gill is based on the 6 macroscopic grades, from 0 (no visual pathology) to 5 (severe visual pathology) (Gill Health Initiative, 2015;Bloecher et al., 2018). It is important to realize that our comparison was done on the groups of tissue that differed only by two grades (PGD 1 and PGD 3) and were in the middle range of the spectrum, with a bias towards less damaged gills. Although our intention was to compare the transcriptome and histopathology of gill tissues with the full spectrum of PGD scores (e.g., PGD 0 and PGD 5) from the same farm location, fish with such a large range in gill gross morphology were not found during the sampling events, not only within the three production sites described in the current study ( Table 2), but also at the four other locations (Scottish Sea Farms) visited by us in years 2017-2019 (R. Bickerdike and S. A. M. Martin, unpublished data). It is unusual for the salmon farmed in the sea to have all 8 gill arches scored as PGD 0 (apart from the first few weeks following the transfer from the freshwater facility), while the macroscopic gill damage classified as PGD 5 is very rare and would not be expected to be found until later in the production cycle, and as a result of a specific acute event or from cumulative significant environmental insults over a period of time (R. Bickerdike, personal communication). The Atlantic salmon genes are grouped by the common human ortholog (BLASTX, E-value < 0.00001, top hit) and sorted from the highest to the lowest Log 2 FC associated with the HGNC identifier. 1 Mean Log 2 FC for A vs C and B vs C comparisons, as detailed in Supplementary Table 7; 2 mean expression of salmon genes mapped to the same human ortholog.
Discrepancies between macroscopic and microscopic examination of the gill tissue have been reported by previous studies, mainly in the context of AGD (Pennacchi et al., 2016). For example, it has been shown that the AGD scoring system (gross morphology) may not be a reliable method of confirming the disease in cases of light severity of AGD (Clark and Nowak, 1999;Zilberg et al., 2001). This is because small AGD-associated lesions, which affect < 10 lamellae and are easy to detect under the microscope, are typically overlooked during the visual inspection of the gill (Adams et al., 2004). As a result, fish exposed to Neoparamoeba perurans and then classified as clear of the AGD symptoms during the macroscopic gill scoring for AGD may in fact need to be re-classified at the level of histopathological examination (Wynne et al., 2008b). Further complexity to the diagnostic problems is added when the gill disease is multifactorial (Wise et al., 2008;Gjessing et al., 2019;Noguera et al., 2019). More broadly, the poor diagnostic and prognostic value of the PGD scores demonstrated in our study is consistent with the limited applicability of the gross morphology to diagnose complex diseases in livestock and humans, which typically need extensive histopathology and molecular profiling to confirm and prognosticate (Hoffmann et al., 2009;Ahmed and Abedalthagafi, 2016;Mobadersany et al., 2018).

Histopathology-Directed Analysis of Gill Transcriptome Between Sites
By sampling sea farmed Atlantic salmon at three production sites (A, B and C) in Scotland, we demonstrated that the gill samples from different sites had different histopathology ( Figure 4B) and different transcriptomic profiles ( Figure 3B), with 240 to 1360 genes identified as differentially expressed between the sites ( Table 4). The drivers of this pronounced site-to-site variability in gill histopathology and transcriptome are unknown and may reflect the differences in (1) fish cohort, including their hatchery origin, past and present husbandry practices and overall health status, (2) local biotic conditions (e.g., diversity of pathogen), (3) local abiotic conditions (e.g., sea water temperature) and (4) interplay and interactions between all the above factors. Identification of these drivers requires more rigorous spatial and temporal sampling regime (Jokinen et al., 2012;Maestrini and Basso, 2018) and is beyond the scope of the current study.
The NMDS plots revealed that the gill samples from site C were different (i.e., clearly separated from sites A and B) not only at the level of their transcriptome (Figure 3B), but also in terms of their histopathology ( Figure 4B). Because the NMDS ordination of the histopathological data is based on the same set of parameters (i.e., LH, LF, LO, CA, in, eg, cc, cd, ib, ch, pp and op) across all fish, lower NMDS1 values of gill samples from site C reflect lower scores of the contributing parameters (lower grade inflammation) than higher NMDS1 values of gill samples from sites A and B (Figure 4B, Supplementary Table 5, and Supplementary Figure 2). In contrast, the NMDS ordination of the transcriptomic data ( Figure 3B) is based on the expression of 500 genes that are most divergent between each pair of fish (amounting to 903 sets of genes for all combinations of 43 fish in total) rather than being performed on the same set of genes across all fish (Robinson et al., 2010). Thus, the higher leading Log 2 FC dim 1 values of gill samples from site C refer to their generally higher magnitude of transcriptomic changes rather than to the specific set of genes that separates them from sites A and B, with latter requiring more detailed analysis of the existing data to perform.
Because gill samples were more different between sites than within sites (Figures 3B, 4B), we performed the histopathology-directed analysis of gill transcriptome on all fish belonging to each site rather than on the individual fish within or between sites (Figure 5). Specifically, we classified all fish from site C as having gills with low histopathology, while all fish from site A and independently all fish from site B as having gills with moderate histopathology. Comparing the two independent sets of the moderately inflamed gill transcriptomes (A and B) with the low inflammation gill transcriptome (C), in the A vs C and B vs C contrasts, led to the identification of the common DEGs that are likely to constitute the gene expression patterns of multifactorial gill inflammation. Similar approaches utilizing independent sampling and site-to-site comparisons to establish common transcriptomic responses to a range of environmental factors have been successfully applied in studies ranging from aquatic toxicology (Oleksiak, 2008;Defo et al., 2018) to human medicine (Mueller et al., 2009;Rohart et al., 2017).

Gene Expression Patterns of Multifactorial Gill Disease
Transcriptome profiling by RNA-seq has become a powerful tool for identification of genes and molecular pathways involved in the progression from health to disease in farmed fish (Sun et al., 2012;Martin and Król, 2017;Houston and Macqueen, 2019;Ronza et al., 2019). Combining RNA-seq assays with an independent sampling regime is especially important for understanding multifactorial diseases, because it allows for the potential identification of the common gene expression patterns of the disease that is essential for developing treatment and prevention strategies (Cookson et al., 2009;Bhuju et al., 2012). This approach has for example been used in the context of non-specific inflammatory diseases in fish intestine (reviewed in Martin et al., 2016). The present study is the first to focus on the transcriptomic patterns of multifactorial gill disease, a condition increasingly prevalent in salmon sea farms.
Transcriptomic comparison of the two independent sets of moderate histopathology gill samples (sites A and B) with the low histopathology gill samples (site C) led to the identification of 462 common DEGs in total, 311 of which were protein-coding transcripts mapped to 235 human orthologs ( Table 5). The identification of these genes as part of the gene expression patterns of multifactorial gill disease was further supported by their uniform direction of changes in the expression levels ( Figure 5). Subsequent functional analysis of these genes by IPA and GO pointed towards a range of immune responses as the most dominant feature of the moderate histopathology gills.
At the gene level, the immune responses were driven by pro-inflammatory cytokines (IL17F, CXCL9, CXCL10, CCL4L1 and TNF superfamily members: LTA and TNFSF14), cytokine receptors (TNFRSF1B, TNFRSF6B and CXCR1) and regulators of cytokine expression and signaling (C1QTNF3, SOCS1 and SOCS3), all of which were upregulated apart from TNFSF14 (Supplementary Table 7). The proteins encoded by IL17 family genes have been shown to stimulate the production of several other cytokines, including IL6 and IL8 in splenocytes of rainbow trout using recombinant IL17A/F2a (Monte et al., 2013) as well as IL1β, IL6, IL8 and TNF-α in head kidney leukocytes of grass carp using recombinant IL17A/F1 (Du et al., 2015). Furthermore, the emerging role of IL17 cytokines in gill mucosal immunity is supported by reports showing differential expression of these genes following the exposure to Aeromonas hydrophila in common carp (Dong et al., 2019) and after challenge with Ichthyophthirius multifiliis in rainbow trout (Syahputra et al., 2019). According to our BLASTX results, CXCL9 represents salmon IL8, whose high expression in our data set (Log 2 FC = 3.4, Supplementary Table 7) is consistent with upregulation of IL17F (Zou and Secombes, 2016). The implication of IL8 in gill immune responses has been previously demonstrated in Atlantic cod (Caipang et al., 2010) and rainbow trout (Olsen et al., 2011;Santana et al., 2016). The protein encoded by LTA (salmon TNF-α2) mediates a large variety of inflammatory, immunostimulatory and antiviral responses in mammals (Upadhyay and Fu, 2013), showing induction at the level of gene expression in gills of Atlantic bluefin tuna during natural infection with Digenea (Pleić et al., 2015), but not in Atlantic salmon following the infection with Neoparamoeba spp. to induce AGD . Cytokines were also overwhelmingly present in our functional analysis of gene expression patterns in the moderately inflamed gills. Specifically, 6 of the top 13 IPA canonical pathways were related to cytokine signaling (Figure 6 and Supplementary Table 8) and 10 of the top 15 IPA upstream regulators were predicted to be cytokines (Figure 7 and Supplementary Table 9). The presence of cytokines among DEGs and also as upstream regulators is consistent with a cytokine signaling cascade, during which one cytokine stimulates its target cells to produce another cytokines (Tisoncik et al., 2012). As a result of cytokine presence among DEGs, 6 of the top 13 IPA downstream effects were predicted to be immunological and inflammatory diseases (Figure 8 and Supplementary Table 10). Altogether, our results identified cytokine-driven immune response as a hallmark of multifactorial gill disease.
Besides the immune response, the multifactorial gill disease from sites A and B had a common transcriptomic profile indicative of tissue damage and repair. Part of the tissue damage could be inflicted by high levels of nitric oxide (NO), suggested by upregulation of NOS2 (Log 2 FC = 2.4, Supplementary Table 7). This gene encodes a NO synthase that is inducible by a combination of lipopolysaccharide (LPS) and certain cytokines (Okamoto et al., 2004). Although NO is produced as a first-line defense against invading pathogens, its strong cytotoxic effects may also damage the tissue of the host (Abramson et al., 2001). Furthermore, high levels of NO may contribute to oxidative stress through production of reactive oxygen species (ROS) (Girouard et al., 2009), which could potentially explain the upregulation of GPX2 in our data set (Log 2 FC = 3.9, Supplementary Table 7). The protein encoded by this gene belongs to the glutathione peroxidase family, members of which catalyze the reduction of organic hydroperoxides and hydrogen peroxide (H 2 O 2 ) by glutathione, and thereby protect cells against oxidative damage (Brigelius-Flohé and Maiorino, 2013). The potential role of GPX2 in the antioxidant defense against increased ROS at gill mucosal surfaces has been recently discussed in the context of exposing Atlantic salmon to peracetic acid (Soleng et al., 2019). It is well established that the repair of inflamed tissue requires efficient removal of damaged cells through controlled cell death (apoptosis) and concurrent cell proliferation to regenerate damaged structures and build up lost tissue, with both processes closely linked to the activity of cysteine-dependent aspartate-directed proteases (caspases) (Fogarty and Bergmann, 2017). In our study, both CASP14 (salmon caspase-14-like) and other genes associated with apoptosis and autophagy (RNF152, CTSL, CTSG, CTSV, RAB32 and TGM2) were upregulated, pointing towards ongoing tissue repair in the moderately inflamed gills (Supplementary Table 7). The increased activity of caspases during gill inflammation has been previously documented in rainbow trout challenged with Aeromonas salmonicida (Rojas et al., 2015) and grass carp after infection with Flavobacterium columnare (Chen et al., 2019). Because of the complexity and size of gill vasculature, the repair of the gill tissue requires extensive vasculogenesis (formation of new blood vessels from vascular precursor cells), angiogenesis (process of outgrowing vessels from the existing vasculature) and arteriogenesis (remodeling of arteries where collateral arterial anastomoses undergo abluminal expansion) (Shi et al., 2017). One of the most potent mediators of new blood vessel formation is angiogenin (Hoang and Raines, 2017), represented in our data set by ANG (salmon angiogenin-1 and ribonuclease-like 3), which was the gene with the highest level of upregulation (Log 2 FC = 4.6, fold change ∼24) in the moderate histopathology vs low histopathology gills ( Table 6 and Supplementary Table 7). This result is consistent with the earlier study conducted by Valdenegro-Vega et al. (2014), showing increased protein levels of angiogenin (fold change ∼12) in the gills of Atlantic salmon following four successive infections with Neoparamoeba perurans. Among other genes involved in angiogenesis were ADM, BMP10 and VCAN, all upregulated in the moderately inflamed gills (Supplementary Table 7). Overall, 5 of top 13 IPA canonical pathways (Figure 6 and Supplementary Table 8) and 5 of top 13 IPA downstream effects (Figure 8 and Supplementary  Table 10) were associated with cell death and proliferation related to tissue damage and repair, highlighting their importance in multifactorial gill disease.

Conclusion and Future Directions
The gill is central for understanding the impacts of climate change on fish health. Recent increases in gill pathologies in sea farmed Atlantic salmon highlight the need to establish the molecular basis of multifactorial gill disease (frequently referred to as PGD, PGI or CGD) to improve diagnosis and preventive measures of this condition. To ensure the multifactorial etiology of gill disease, we sampled Atlantic salmon from three different production sites in Scotland and then examined the gill tissue at three different levels of organization: gross morphology with the use of PGD scores (macroscopic examination), whole transcriptome (gene expression by RNA-seq) and histopathology (microscopic examination). By exploring the association between gill transcriptome and gill gross morphology (PGD scores), we found that the PGD scores were less effective in providing a graded assessment of gill health status than expected, and they did not convey any information about the underlying pathology and/or tissue deterioration. In contrast, integration of the gill RNA-seq data with the gill histopathology enabled us to identify common gene expression patterns associated with multifactorial gill disease. We demonstrated that the gene expression patterns associated with multifactorial gill disease were dominated by two processes: a range of immune responses driven by pro-inflammatory cytokines and the events associated with tissue damage and repair, driven by caspases and angiogenin.
Previous studies of the gill inflammation have typically focused on single-cause pathologies (such as AGD), in the experiments performed on fish exposed to the controlled environment of the closed-system tanks. Although these studies are very important from a mechanistic point of view of specific pathologies, they may have limited relevance for the fish facing a changing environment at farm locations, e.g., rising sea water temperature, deoxygenation and surges of existing and emerging pathogens with the complexity of co-infections, coupled with husbandry practices. Performing more transcriptomic studies in the field rather than in the lab would benefit both academia and industry. The multisite approach is also important because in this study, we demonstrated large and significant effects of sampling site (sensu lato) on gill transcriptome and histopathology. The drivers of this site-to-site variability are presently unknown and require more specific sampling regime.
Our histopathology-directed analysis of gill transcriptomes identified in total 462 genes that we claim to constitute the gene expression profile of multifactorial gill inflammation, including 354 protein-coding genes, 35 immunoglobulin gene segments, 15 pseudogenes and 58 non-coding RNAs. It is, however, important to remember that our analysis was based on the gill samples from three production sites in one specific part of the world (coastal waters of Scotland), with sampling events covering October, November and March. Substantial amount of work is therefore needed to test the association of these genes with multifactorial gill diseases at different times of year in Scotland and worldwide. The diagnostic and therapeutic value of these transcripts is currently unknown and require further studies.

DATA AVAILABILITY STATEMENT
The RNA-seq data have been deposited in the ArrayExpress repository (http://www.ebi.ac.uk/arrayexpress/) under accession number E-MTAB-8855. The R and Python scripts are available on request to AD and SS, respectively.

ETHICS STATEMENT
Ethical review and approval was not required for the animal study because the sampling of fish was carried out under established protocols for routine health assessments in accordance with RSPCA Assured Welfare Standards for farmed Atlantic salmon and under supervision of the company veterinarian.

AUTHOR CONTRIBUTIONS
SM, AD, EK, VV, KG, and RB conceived and designed the study, interpreted results, and gave final approval of the manuscript. All authors (apart from SS and AD) were involved in sampling. EK extracted RNA, performed IPA analysis, and drafted the manuscript with AD, who oversaw bioinformatics and statistical analysis. SS mapped salmon genes to human orthologs to enable functional analysis of gene expression. PN performed histopathological examination of samples, generated scores, and interpreted results. EC provided support in data exploration. All authors discussed and commented on the manuscript.

FUNDING
The study was supported by the Scottish Aquaculture Innovation Centre (SAIC grant SL 2017 08, 'Nutritional Aspects of Gill Disease in Atlantic Salmon'). The authors declare that this study received funding from BioMar AS and Scottish Sea Farms (SSF). The funders had the following involvement with the study: contribution to the experimental design, selection of sampling sites, participation in sampling, contribution to the interpretation of the results, and approval of the manuscript.

ACKNOWLEDGMENTS
We thank SSF farm personnel for accommodating our research, performing gross morphology scoring and helping with sampling. Edinburgh Genomics is partly supported through core grants from NERC (R8/H10/56), MRC (MR/K001744/1), and BBSRC (BB/J004243/1). We also thank two reviewers for their comments on the earlier draft of the article. EK dedicates this paper to her late mother, Irena Król.